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ABSTRACT 

The outward migration of a pair of resonant-orbit planets, driven by tidal interactions with a gas- 
dominated disk, is studied in the context of evolved solar nebula models. The planets' masses, Mi 
and M%, correspond to those of Jupiter and Saturn. Hydrodynamical calculations in two and three 
dimensions are used to quantify the migration rates and analyze the conditions under which the 
outward migration mechanism may operate. The planets are taken to be fully formed after 10 6 and 
before 3 x 10 6 years. The orbital evolution of the planets in an evolving disk is then calculated 
until the disk's gas is completely dissipated. Orbital locking in the 3:2 mean motion resonance may 
lead to outward migration under appropriate conditions of disk viscosity and temperature. However, 
resonance locking does not necessarily result in outward migration. This is the case, for example, 
if convergent migration leads to locking in the 2:1 mean motion resonance, as post-formation disk 
conditions seem to suggest. Accretion of gas on the planets may deactivate the outward migration 
mechanism by raising the mass ratio M2/M1 and/or by reducing the accretion rate toward the star, 
hence depleting the inner disk. For migrating planets locked in the 3:2 mean motion resonance, there 
are stalling radii that depend on disk viscosity and on stellar irradiation, when it determines the disk's 
thermal balance. Planets locked in the 3:2 orbital resonance that start moving outward from within 
1-2 AU may reach beyond ss 5 AU only under favorable conditions. However, within the explored 
space of disk parameters, only a small fraction - less than a few percent - of the models predict that 
the interior planet reaches beyond ps 4 AU. 

Keywords: accretion, accretion disks — hydrodynamics — methods: numerical — planet-disk inter- 
actions — planets and satellites: formation — protoplanetary disks 



1. INTRODUCTION 

The architecture of the solar system bears some evi- 
dence that Jupiter and Saturn may have been closer to 
each other in the past (e.g., Malhotra 1993, 1995; Tsi- 
ganis et al. 2005; Morbidclli ct al. 2005; Gomes et al. 
2005) . They later moved away from each other because of 
gravitational interactions with the remnants of the disk 
of planetesimals from which these planets had formed 
(e.g., Fernandez & Ip 1984; Hahn & Malhotra 1999). The 
planetesimal-driven migration of Jupiter and Saturn oc- 
curred relatively late, after the gaseous component of the 
solar nebula had dispersed, and the extent of their ra- 
dial displacements was probably less than ~ 1 AU (e.g., 
Franklin et al. 2004; Minton & Malhotra 2009). 

Recently, Walsh et al. (2011) proposed a scenario in 
which orbital migration of Jupiter and Saturn occurred 
much earlier in the solar system history and was driven 
by tidal torques in a gas-dominated nebula. The progen- 
itors of Jupiter and Saturn underwent rapid convergent 
migration toward the Sun, until Saturn became trapped 
in the 2:3 mean motion resonance with Jupiter. By 
that time and under the applied conditions, Jupiter had 
reached ps 1.5 AU and Saturn ps 2.0 AU. Once the reso- 
nant configuration was established, the planets reversed 
the direction of motion and began migrating outward, 
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preserving the 2:3 commcnsurability. This scenario may 
help explain some features of the inner solar system, in- 
cluding the Mars-to-Earth mass ratio and the radial vari- 
ation of composition in the asteroid belt (see Walsh et al. 
2011, for details). 

The outward migration is a direct result of the "com- 
pact" orbital configuration. Qualitatively, the negative 
torque balance that would result for a single-planet is 
tipped in favor of the positive torque (from the inner 
disk) because the negative torque (from the outer disk) 
is abated by a local reduction of the surface density. This 
situation requires that the planets be massive enough to 
significantly perturb, via tidal interaction, the disk's sur- 
face density and that their density gaps overlap. These 
requirements are typically realized if the orbital separa- 
tion is at most several times the sum of the planets' Hill 
radii. Therefore, depending on the masses, a (near) 3:2 
commensurability is favorable to sustain outward migra- 
tion of a Jupiter-Saturn pair, whereas for more massive 
planets, by a factor of about three, a (near) 2:1 commen- 
surability may promote outward migration. 

A study by Pierens & Raymond (2011) lends support, 
under appropriate conditions, to the inward-outward mi- 
gration scenario of the Jupiter-Saturn system proposed 
by Walsh et al. (2011). One scope of this paper is to 
revisit this idea in the context of evolved models of a 
gas-dominated solar nebula. In particular, we concen- 
trate on the outward migration of a pair of giant planets, 
whose masses correspond to those of Jupiter and Saturn, 
after their orbits become locked in the 3:2 mean motion 
resonance, compatibly with the formation timescales of 
both Jupiter and Saturn, estimated from core-nucleated 
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accretion models. 

We also wish to provide some constraints on the range 
of radial migration of Jupiter (and Saturn), as a func- 
tion of the solar nebula properties, under the assumption 
that the 3:2 orbital resonance is maintained throughout 
the disk's evolution. Conditions that may break the res- 
onance locking between the two planets or that may in- 
hibit or prevent outward migration are analyzed as well. 
In particular, we focus on the process of gas accretion 
that, on one hand, may alter the planets' mass ratio and, 
on the other, may reduce the disk density inside the or- 
bit of the interior planet. Both effects act to change the 
balance of the torques exerted on the planets. In addi- 
tion, we examine the disk conditions under which con- 
vergent migration leads to capture of the exterior planet 
in the 1:2 orbital resonance with the interior planet, a 
configuration that does not promote outward migration 
of a Jupiter-Saturn pair, and which may leave the plan- 
ets stranded in the inner disk region. The possibility 
that Saturn forms within the 1:2 commensurability with 
Jupiter is also analyzed. 

The layout of the paper is as follows. In Section 2, we 
describe dynamics and thermodynamics of disk models 
and report on their evolution. In Section 3, the tidal 
interaction calculations in two and three dimensions are 
presented, along with the calculations of the migration 
rates of a 3:2 resonant-orbit pair. Section 4 is dedicated 
to the long-term orbital evolution of two planets locked 
in the 3:2 orbital resonance. Two possible effects of gas 
accretion are analyzed in Sections 5 and 6, while condi- 
tions for capture in the 2:1 mean motion resonance and 
some related issues are examined in Section 7. Section 8 
contains the discussion and the summary of the results. 

2. LONG-TERM DISK EVOLUTION MODELS 

In this section we describe the dynamics and thermo- 
dynamics of solar nebula models. For tested parameters, 
we report on the disk evolution until the gas is almost 
entirely dispersed, that is until the disk mass, Md, is less 
than 10 -5 times the mass of the star. By assumption, 
successful sets of parameters representing a solar neb- 
ula model are those that provide a disk lifetime, td, no 
greater than ~ 2 x 10 7 years. Although the gas mass of 
disks is notoriously difficult to ascertain, according to ob- 
servations (see, e.g., reviews by Roberge & Kamp 2010; 
Williams & Cieza 2011), the presence of gas in the in- 
ner regions of protoplanetary disks appears to last < 10 7 
years (see also Haisch et al. 2001). 

2.1. Disk Dynamics 

Consider a gaseous disk orbiting a central star of mass 
M s . In the framework of one-dimensional (ID) mod- 
eling, we assume azimuthal symmetry around the star 
and use vertically averaged quantities as a function of 
the radial distance r. For the current purposes, we as- 
sume that the evolution of the disk is driven by viscous 
torques, l~ v , and wind dispersal, M w at the disk's sur- 
face. The torque exerted on a disk ring of radius r, by 
material orbiting inside the ring, is 71/ = —2nr 3 vY.dfl/dr 
(Lyndcn-Bcll & Pringlc 1974), where v is the kinematic 
viscosity of the gas, £ the surface density, and the an- 
gular velocity. If f2 is identified as the Keplerian velocity 
(i.e., if effects of gas and magnetic pressure gradients arc 



neglected), then T v — 'invTM, where T~L = r 2 £l is the spe- 
cific angular momentum of the gas. Along with viscous 
diffusion, the disk is dispersed by a wind, whose origin 
is gas photo-evaporation from the disk surface produced 
by photons emitted by the central star. Hence, we write 
M w = 2ir J S po r(ir, where S pe is the mass per unit sur- 
face area and unit time removed from the disk. 
The continuity equation for the disk requires that 
19 



(1) 



where u r is radial velocity of the gas. On the left-hand 
side, one can recognize the mass per unit time flowing 
through a circumference of radius J- = 2irrY,u r . By 
using the relation T = —dTv/dH (see Lyndon-Bell & 
Pringle 1974) and since we assume Keplerian rotation 
{dH/dr = rfl/2), Equation (1) becomes 



^_(S + S pe )--I- — 



0. 



(2) 



To seek for numerical solutions of Equation (2), it is con- 
venient to use Ti as independent variable and S = "H 3 £ 
as dependent variable and then solve 



S ve )-\{GM S ?-^ 
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In the above equation, G is the gravitational constant 
and S pc = - H 3 S po . Note that quantity S /Ti 2 is the an- 
gular momentum per unit surface area. In writing Equa- 
tion (1), we neglected the effects of the star's growth, 
which would introduce a term on the right-hand side of 
order EM S /M S (see Ruden & Pollack 1991). Given the 
initial values of M s /M s considered here (see Section 2.3) 
and the decline of M s with time, this term would affect 
£ only over a time scale of order 10 years, or longer. 

Photo-evaporation involves contributions from far- 
ultraviolet (FUV), extreme- ultraviolet (EUV), and X-ray 
radiation emitted by the star (see Dullemond et al. 2007; 
Clarke 2011, and references therein). FUV radiation may 
be especially important in removing gas at large dis- 
tances from the star, reducing the gas supply to the inner 
parts of the disk. However, a self-consistent calculation of 
FUV photo-evaporation rates requires solving for the de- 
tailed vertical structure of the disk (e.g., Gorti et al. 2009; 
Gorti & Hollenbach 2009). Photo-evaporation by EUV 
photons is more tractable since they ionize hydrogen at 
the very upper layers of the disk. Here we follow a simple 
approach and adopt the formulation of the EUV photo- 
evaporation rate proposed by Dullemond et al. (2007): 



(^) 5/2 



for r > r„. 



(4) 



The radius r g w 10 (M s /Mq) AU is the gravitational 
radius, beyond which gas at the disk surface is unbound 
(see, e.g., Armitage 2011, and references therein). The 
photo-evaporation rate at r g is 



E!L = 1.16 x 10 



-11 




(5) 



where fn is the rate of EUV ionizing photons emitted 
by the star in units of 10 41 s . The total mass loss rate 
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due to photo-evaporation is found by integrating Equa- 
tion (4) over the entire disk according to the definition 
given above, hence M w = (0.55977\/e + 47r) r 2 E| e or 



M w = 1.56 x 10" 10 W/ 4 i 



1AU 



Meyr" 1 . 



(6) 



The maximum of E pe occurs at r 



r g /4. Locally, 

gas is removed via photo-evaporation and supplied by 
viscous diffusion, i.e., accretion through the disk, M = 
—T (note that T is positive for an outward transfer of 
mass). Recalling the relations reported above, we can 
write dTv/d'H = STtd{vYfH) / 'dH, and thus 



M = 3tt 



2r|>S) 
or 



(7) 



For vYi nearly independent of r, i.e., in a stationary disk 
(Pringle 1981, see also Equation 1 with the right-hand 
side set to zero), M = 37wE is nearly constant through- 
out the disk. Therefore, if M s = 1M Q , we expect gas 
depletion induced by photo-evaporation to occur first 
around ~ 3 AU. 

2.2. Disk Thermodynamics 

In order to determine the thermal energy budget of 
the disk during its evolution, we assume that there is a 
balance among three terms: viscous heating, irradiation 
heating by the central star, and radiative cooling from 
the disk's surface. Viscous dissipation produces an en- 
ergy flux equal to Q u = uY, (rdtt/dr) (see, e.g., Mihalas 
& Weibel Mihalas 1999), which in case of the Keplerian 
rotation becomes 



Q u = ^Eft 2 . 
4 



(8) 

Since Q u oc 1/r 3 , for a disk with d{uT,)/dr w 0, vis- 
cous dissipation becomes an ever less important source 
of heating as the distance from the star increases. 

We follow the formulation of Hubcny (1990) for an ir- 
radiated disk and write the energy flux escaping from 
both sides of the disk surface as 



'3 1 1 

Qcool = 2(7SB T [ -Tr + - + -j— - 



(9) 



whereas the heating flux arising from stellar irradiation 
can be written as 



Q irr — 2<jsb T irl 



TR 



4r P y 



(10) 



In the above equations, ctsb is the Stcfan-Boltzmann con- 
stant, T the mid-plane temperature, and T lrr the irradi- 
ation temperature. Note that, for an irradiated disk, the 
constant in parenthesis on the right-hand side of Equa- 
tion (9) is generally slightly different from that of a non- 
irradiated disk (compare with Equation 14 of D'Angelo 
et al. 2003). As in Menou & Goodman (2004), we set 



T 



(l-e)T s 4 



Wq, 



(11) 



where e is a measure of the disk's albedo, for which we 
adopt the value 1/2, and T s and R s are the effective 
temperature and radius of the star, respectively. This 



interpretation of the irradiation temperature, however, 
neglects the contribution of luminosity released by stellar 
accretion (e.g., Hartmann et al. 2011). In an actively 
accreting disk, quantity T s 4 should be replaced with — 
+ T^cci where T^ cc quantifies the luminosity due to 
accretion L acc = GM S M S /(2R S ) (Pringle 1981) and thus 



rjil -*- 




(12) 



where the accretion rate M s , computed as —d%/dH (sec 
Section 2.1) at the disk's inner radius, varies with time. 

The quantity Wq in Equation (11) is a geometrical fac- 
tor that accounts for illumination of disk portions close 
to (first term) and far from (second term) the star (see 
Chiang & Goldreich 1997) 



Wq = 0.4 



R, 



H 

r 



dbiH 
G?lnr 



(13) 



T he adiabatic scale-height of the disk, H = 
7 UbT / (/xttt-h) /fi, is derived from the requirement of 
vertical hydrostatic equilibrium. The adiabatic index, 7, 
is 1.4, the mean molecular wight, fj,, is 2.39, fee is the 
Boltzmann constant, and mn the hydrogen mass. 

If the second term on the right-had side of Equa- 
tion (13) is negative, the disk is self-shadowed and that 
term should be dropped. A self-consistent calculation 
of this term from ID, vertically averaged models may 
lead to numerical instabilities (see, e.g., Hueso & Guil- 
lot 2005) . In fact, meaningful determinations of this 
term involve solving for the vertical thermal structure of 
the disk. Therefore, the last term in parenthesis on the 
right-hand side of Equation (13) is written as rj and ap- 
proximated to 2/7 (see, e.g., D'Alessio et al. 1998; Menou 
& Goodman 2004; Hueso k Guillot 2005; Rafikov & De 
Colle 2006). 

The optical depths tr = krE/2 and rp = KpE/2 in 
Equations (9) and (10) are based, respectively, on Rosse- 
land (kr.) and Planck (kp) mean opacities. Both kr 
and Kp depend on T and the mass density p = T,/(2H). 
We adopt grain opacities from Pollack et al. (1994), at 
temperatures below the vaporization temperatures of sili- 
cates, and gas opacities from Ferguson et al. (2005) for so- 
lar abundances, when all grain species have evaporated. 

The thermal energy budget is given by 



= 0. 



(14) 



Note that if Q u <C Qi rr , a situation that may occur in an 
evolved disk, Equation (14) results in a gas temperature 
T = Tin-, that is 



T 




e)W c 



11/4 



(15) 



The factor Wq is typically a weakly dependent function 
of T. If Wq is a constant, then T cx r" 1 / 2 . If Wq oc R s /r 
(e.g., at radii r ~ R s ), then T cx r~V 4 H Wq oc H/r 
(as we assume for r ^> R s ), then W G oc T 1 / 8 , the 
temperature is T oc r -3 / 7 (sec also Chambers 2009), and 
the disk's aspect ratio is 



HV = (jk^Yf R 

r J \ ^m H 



>{^ 2 - (16) 
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The choice of the parameter r\ may have some impact 
on the disk's thermal budget, yet Equation (16) suggests 
that this impact is low. 

2.3. Numerical Procedures and Parameters 

Equation (3) is evolved in time using an implicit 
numerical scheme, which avoids the sometimes pro- 
hibitively short time steps required by an explicit ap- 
proach, especially when the inner disk radius extends 
very close to the star (see Bath & Pringlc 1981). We 
either use a second-order Crank-Nicolson method (e.g., 
Press et al. 1992) or a fourth/fifth-order Dormand- 
Prince method with an adaptive step-size control based 
on the global accuracy of the solution (Haircr ct al. 1993). 
In the latter case, the evaluation of derivatives (in the 
Runge-Kutta sequence) is performed by means of a back- 
ward Euler (implicit) method. A zero-torque boundary 
condition, 5 = 0, is applied at the disk's inner edge. At 
the outer edge, the applied boundary condition is such 
that dM/dH = d 2 %/dH 2 is constant. Figure 1 (left) 
shows a comparison between numerical (lines) and ana- 
lytic (circles) solutions of Equation (3) (see figure's cap- 
tion for details). 

Equation (14) is solved for the mid-plane temperature, 
T, at each radius, using a root-finding algorithm based 
on the Brent's method (Brent 1973). Convergence of 
the root-finding process is achieved within a tolerance of 
10 -3 K. An iterative procedure is implemented for each 
determination of T, so that the applied value of H and 
that corresponding to the converged temperature do not 
differ by more than 1%. In Figure 1 (right), the evolu- 
tion of temperature is shown for the disk considered in 
the left panel. In this test, we set Wg = 0.05, so that 
temperature evolves toward that in Equation (15), indi- 
cated as filled circles. The temperature profiles show ma- 
jor opacity transitions at T w 160, 420, 680, and 1400 K, 
caused by vaporization of, respectively, water ice, refrac- 
tory organics, troilite, and silicate grains (see Pollack 
et al. 1994). Note that heating via viscous dissipation 
is basically confined within ~ 10 AU (see Section 2.2). 

The solar nebula extends from r m = 0.01 AU to r out = 
1850 AU and is discretized over 10000 grid points. The 
large outer radius is chosen to not interfere with viscous 
spreading of the disk. The numerical resolution is vari- 
able and such that Ar/r ~ 1.2 x 10~ 3 . In this study, we 
assume that M s = 1M Q , T s = 4280 K, and R s = 2R Q 
(Siess et al. 2000). The initial surface density distribu- 
tion of the gas obeys the relation S = E^ri/r)' 3 , where 
13 = 1/2, 1, or 3/2, within at least ~ 10AU. Farther 
away from the star, £ is exponentially tapered. The ex- 
tremes of the "slope" (3 bracket values derived for the so- 
lar nebula by Davis (2005) and by Weidcnschilling (1977) 
and Hayashi (1981). The quantity the surface den- 
sity at ri = 1 AU, is such that the initial disk mass is 
Mg ~ 0.022, 0.044, or 0.088 M s . (These will be re- 
garded as nominal values. The total initial disk mass 
differs somewhat for the different values of (3 because 
of the tapering procedure). The photo-evaporation rate 
(Equation 4) is specified by imposing fu in Equation (5). 
Here we use / 4i = 0.1, 1, 10, 100, and 1000 (Alexander 
et al. 2005). The kinematic viscosity is v = v\(r/r\f and 
v x = 4x 10~ 6 , 8x 10" 6 , and 1.6 x 10" 5 rf Qi, where fii is 
the rotation rate at r = r\. As a reference, in a disk with 
constant aspect ratio H/r = 0.04, i/i = 8x r\£l\ 



corresponds to a turbulence parameter (Shakura & Syun- 
yaev 1973) at — 0.005. The initial accretion rate onto 
the star ranges from a few times 10 -8 to a few times 
10 -7 Afoyr -1 . For comparison, the mass loss rate in 
Equation (6) is between ~ 10" 10 and ~ 10~ 8 M yr _1 

2.4. Model Results 

The majority of disk models have an initial gas inven- 
tory of at least ~ 0.02 M Q within a distance of 40 AU 
from the Sun, as required by a canonical minimum mass 
solar nebula (MMSN; e.g., Weidenschilling 1977; Hayashi 
1981). This value is also consistent with the more recent 
MMSN model adopted by Chiang & Youdin (2010). Due 
to the steepness of the surface density, disk models with 
the lowest initial mass and parameter ,8 = 3/2 have only 
0.01 M Q worth of gas within 40 AU of the Sun. 

Gas is removed via the combined action of accretion 
onto the star, M (Equation 7), and photo-evaporation, 
M w (Equation 6). In particular, Equation (6) sets an 
upper limit to dispersal timescale, ttj, ranging from 

- 1.4 Myr for MS ~ 0.022 M s (when f 41 = 1000) to 

- 560 Myr for Mg ~ 0.088 M s (when / 41 = 0.1). For 
computational purposes, td is defined as the time past 
which M u < 10~ 5 M s . 

The evolution of the disk mass for some selected cases 
is illustrated in Figure 2 for each reference viscosity (see 
figure's caption for details). A complete list of the disk 
lifetimes is reported in Table 1. The behavior of the 
disk mass as a function of time, for the different sur- 
face densities, can be qualitatively understood in terms 
of viscous evolution by means of the analytic solutions 
of Lynden-Bell & Pringle (1974, their Section 3.3): for 
equally massive disks, the more compact the disk is (i.e., 
the larger /?), the more rapidly Mb reduces initially. By 
a somewhat conservative assumption, as discussed above, 
disks that survive beyond 20 Myr are discarded and will 
not be given any further consideration. This is the case, 
for example, for all models with a photo-ionizing rate 
characterized by fn < 1 and the flattest initial surface 
density (j3 = 1/2). Models of disks surviving less than 
1 Myr will also be discarded based on considerations on 
planet formation timcscales, as explained in Section 4. 

A quantity of primary importance for planetary migra- 
tion is the average surface density around the planet's or- 
bit. In Figure 3 (left panels), the evolution of £ is shown 
for cases with different values of parameters (3 and /41 . 
As anticipated at the end of Section 2.1, once the accre- 
tion rate drops below some threshold, photo-evaporation 
produces a gap in the surface density at a radial distance 
of a few AU. Then the disk inside the gap is removed 
by viscous diffusion on a timescale of order r 2 /(2ttv) or- 
bital periods. We shall see in the Sections 3.2 and 3.3 
that the disk's aspect ratio has also a large impact on 
the rates and direction of migration. In the right pan- 
els of Figure 3, H/r is plotted at reference times for the 
same models as in the left panels. Once E becomes small 
enough and the viscous heating term Q v in Equation (14) 
becomes unimportant, the disk temperature, and hence 
H/r (see Section 2.2), is dictated only by the stellar ir- 
radiation temperature (Equation 11). 

The surface density at 1 AU, £1, vs. time is illustrated 
in Figure 4 for selected models from Table 1. Thin and 
thick lines refer to smallest and largest values of v\, re- 
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Figure 1. Left: evolution of a viscous disk obtained by solving Equation (3), with iS po = 0, by means of the Dormand-Prince method. 
The initial condition (open circles) is the analytic solution of Lynden-Bell & Pringle (1974) for a disk with no central couple, Mj3 = 0.1 Mq, 
M a = 10 — 7 Mq yr —1 , and u = 8 X 10 — 6 rf n Q- ln . The solid lines represent the numerical solution at different times and the filled circles are 
computed using the analytic solution. Right: evolution of temperature obtained by solving Equation (14) for the disk in the left panel. 
Times in the legend are in years. For testing purposes, Wg is set equal to 0.05. The temperature predicted by Equation (15) is indicated 
by filled circles. 
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Lifetimes from Disk Models 
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a Time past which M D < 10~ 5 M s , in units of Myr. 
b Initial "slope" of the disk's surface density. 

c Kinematic viscosity at r\ = 1 AU in units of r\ Q\ = (G M s ri) 1/2 . 

d Rate of EUV ionizing photons emitted by the star in units of 10 41 s — 1 (see Equation 5). 



spectively. To obtain a better statistical characterization 
of the disk density and temperature, for each pair of pa- 
rameters (M D , (3) listed in Table 1, v\ and /41 are varied 
randomly in the corresponding ranges indicated in the 
table, for a total of more than 2000 realizations. The 
histogram in Figure 5 (left) shows that after ~ 1 Myr 
the value of Si is 150 gem -2 , or less, in ~ 80% of the 
models that may represent the solar nebula. The right 
panel of Figure 5 illustrates the occurrence frequency of 
the mid-plane temperature at 1 AU. As a reference, the 
equilibrium temperature established by stellar irradia- 
tion alone, i.e., Equation (15), is T\ rts 100 K. 

3. TIDAL INTERACTIONS OF JUPITER AND 
SATURN WITH THE DISK 

The evolution of the thermodynamical quantities 
(principally E, T, and H/r) of solar nebula models, ob- 
tained in the previous section, can be used to evaluate 
the range of orbital migration of a pair of planets, over 
the disk lifetime, once appropriate migration rates are 
supplied. In this section we derive such rates. 



3.1. 2D and 3D Hydrodynamical Calculations 

The migration of a Jupiter-Saturn pair in a gaseous 
disk is evaluated by using a combination of 2D and 
3D hydrodynamical calculations of tidal interactions be- 
tween the planets and the disk. We adopt a reference 
frame {0; r, 0, </>} with origin, O, fixed on the star, radius 
ranging from 0.25 to 7 AU, and azimuth varying between 
and 2tt. In the 2D disk approximation, the co-latitude 
angle 9 is equal to 7r/2, whereas it varies from 9 mm to 
7r/2 in a 3D disk. In the latter case, the disk opening 
angle, O m ini is such that the disk's vertical extent locally 
comprises at least three pressure scale-heights, H. Mir- 
ror symmetry with respect to the 6 = ir/2 plane is im- 
posed on account of the planets orbiting in this plane of 
symmetry. The surface density E has initially a dimen- 
sionless gradient dlnE/rflnr = —1/2. We work in the 
assumption that the disk is locally isothermal, i.e., the 
temperature depends only on r, and that H/r is a con- 
stant. The gas pressure is therefore proportional either 
to E/r (2D) or to p/r (3D), where p is the mass density. 
It is further assumed that the kinematic viscosity of the 
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Figure 2. Mass evolution for disks with initial (nominal) masses M° ~ 0.022 M s (left) and 0.088 M a (right). The initial S has /3 = 1/2 



(top), 1 (center), and 3/2 (bottom). Thin and thick lines represent models with viscosity i/i=4x 10 — 6 and 1.6 X 10 — 5 
Different line colors correspond to different rates of EUV ionizing photons emitted by the star in units of 10 41 s — 1 
legends. 



disk, v, is constant throughout the disk. 

The coordinate system rotates about the axis perpen- 
dicular to the planets' orbital plane (0 = axis) at a 
variable rotation speed, fi/ = H/(£). Both Clt and Clf 
are imposed by the requirement that the (relative) az- 
imuthal position of the interior planet, </>i, remains, con- 
stant in time and the (relative) angular velocity, </>i, is 
zero (for details about the procedure, see D'Angelo et al. 
2005). 

Naming ri and r 2 the vector positions of the interior 
(Jupiter) and exterior planet (Saturn), respectively, the 
gravitational potential in the disk is 



GAh 



-r ■ n 



GM 2 



-r ■ r 2 , 



' Hi , respectively, 
as reported in the 



(17) 



GM 2 



y/\T- ri \ 2 + e'i v/|r-r 2 | 2 +e 



which accounts for the contributions of non-inertial terms 
due to the reference frame being centered on the star (see 
Nelson et al. 2000). The potential softening lengths E\ 
and e 2 are set equal to 1/4 (or 1/7, in some calculations) 
times the Hill radius, i?H, of the corresponding planet. It 
is worth stressing that the argument according to which 
e should be a fraction of the disk's scale-height, H , in 
the 2D geometry (e.g., Masset 2002; MMer et al. 2012) 
applies to fully embedded planets, when i?n < H. If 
i?H ft H , as in all our 2D calculations, the local disk 
scale-height depends also on the gravity of the planet 
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Figure 3. Evolution of the disk's surface density (left) and disk's relative thickness (right) of disk models with v\ = 
MK ~ 0.022 M 3 , / 4 i = 10, and /3 = 1/2, (top), v\ = 8xW~ 6 rf Hi, Mg ~ 0.088 M s , ,f 4 i = 100, and /3 = 1 (center), and vy 
MS ~ 0.044 M s , fa = 1, and $ = 3/2 (bottom). Times indicated in the legend arc in Myr. 



S X 10" 
: 4x10" 



6 rf f2i, 



itself (e.g., D'Angelo et al. 2003). In such cases, one 
physical constrain on e is that it should be smaller than 
the radius over which gas is effectively bound to (i.e., it 
rotates about) the planet (~ -Rh/3, see, e.g., D'Angelo 
et al. 2003). 

The Navier-Stokes equations that characterize the disk 
evolution are solved by means of the finite-difference 
code described in D'Angelo et al. (2005, and references 
therein) with modifications detailed below. The disk is 
discretized in 678 x 16 x 700 grid zones, in r, 6, and 
4>, respectively (and 678 x 700 in 2D). Calculations were 
also performed at a higher resolution of 1353 x 28 x 2096. 
Comparisons of the evolution of the planets' orbital el- 
ements at these two resolutions yield good agreement. 
We apply the wave-damping boundary conditions of dc 



Val-Borro et al. (2006) within r = 0.3 and beyond r = 
6.65 AU, which are appropriate for planets far enough 
from the boundaries. 

We implemented an orbital advection algorithm along 
the lines of the FARGO algorithm of Masset (2000) (see 
also Kley et al. 2009). These types of algorithms exploit 
periodicity properties of the flow, as those naturally oc- 
curring in the azimuthal direction of a disk. (Note that 
these algorithms can also be applied to local disk simula- 
tions, if periodicity is imposed at the patch boundaries, 
see Gammie 2001). As demonstrated by Masset (2000), 
when the highest velocity component is along the peri- 
odic direction, in a 2D disk such techniques can increase 
the time step limit required by the Courant-Friedrichs- 
Lewy condition (see, e.g., Stone & Norman 1992), rela- 
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Figure 4. Evolution of the disk's surface density at 1 AU. Thin and thick lines represent models with viscosity v\ = 4 X 10 and 
1.6 X 10 — 5 r\ Oi, respectively. Different line colors correspond to different disk's initial masses. The value of the EUV ionizing photon rate, 
/41, and of the initial surface density gradient, /3, are given in each panel. 



tive to a standard advection scheme, by factors ~ 10, 
or larger. In a 3D disk, the gain may critically de- 
pend on the numerical resolution in the vertical (9) direc- 
tion 5 . The implementation requires care when handling 
the transport of quantities defined on staggered meshes. 
Kley et al. (2009) use split cells, apply the transport al- 
gorithm to each part, and then recombine the partial in- 
formation to reconstruct the full transport. Unlike them, 
we define the auxiliary variables required in the procc- 

5 Since a disk is typically thin, very high resolutions can be 
more easily achieved in the vertical direction. In a viscous disk, 
the time step constraint required by the diffusion part of Navier- 
Stokes equations scales as the square of the gird spacing. This 
requirement, at high resolution and high viscosity, can severely 
reduce (or even nullify) the benefits of orbital advection. 



dure (see Masset 2000, for details) on the same staggered 
meshes as the transport quantities are defined, wherever 
they are necessary. This approach requires more copies of 
the standard auxiliary variables to be defined, and hence 
more memory storage, but offers the advantage that the 
advection of all hydrodynamical quantities can be per- 
formed in a single step. Contrary to the implementations 
of both Masset (2000) and Kley ct al. (2009), the algo- 
rithm used here avoids any directional bias, maintaining 
the full symmetry of the advection scheme, by using a 
sequence that alternates the transport among directions 
(see Stone & Norman 1992). 

The equations of motion of two planets orbiting in a 
disk around a star, written in a reference frame rotating 
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Figure 5. Histograms of the surface density (left) and mid-planc temperature (right) at 1 AU, at time t ss 1 Myr for disk models whose 
lifetime Trj is longer than 1 and shorter than about 20 Myr. The histograms are based on over 2000 models using a random selection of v\ 
and /41, for each pair (M^,/3) listed in Table 1, in the respective ranges indicated in the table. 



at variable angular speed, are 

G(M S +Mi) GM 2 GM 2 
i"i = ^ ri r 2 3— r i2 +Ai-A s 

— flfx(flfxri) — 2 flfxri — flj xri 



1*2 : 



G(M S +M 2 ) 



-r 2 



GM\ 



-it 



GM X 



ri2 + A 2 



'2 ' 1 '12 

-f2^x(n^xr 2 ) — 2 fiy xf 2 — fi/xr 2 , 



(18) 

-As 
(19) 

where ir 2 = it — r 2 . Note that, since the origin of the 
coordinate system is on the star, Equations (18) and (19) 
include the forces per unit mass exerted on the star by the 
planets. The gravitational acceleration terms imposed by 
the disk, Ai, A 2 , and A s , are defined by Equations (8) 
and (9) of D'Angelo et al. (2005) and updated every hy- 
drodynamical time step, At. Equations (18) and (19) are 
integrated numerically over At by means of a high-order 
Gragg-Bulirsch-Stoer extrapolation algorithm with order 
and step-size control (Hairer et al. 1993). The algorithm 
chooses automatically a suitable order at each (internal) 
step, which basically depends on the required tolerance 
of the solution error. We set a relative tolerance of 10 ~ 9 
and an absolute tolerance of 10~ 14 . 

3.2. Torque Calculations and Outward Migration 

The basic mechanism that may allow a pair of 
resonant-orbit planets to experience a positive torque ex- 
erted by a gaseous disk and migrate outward was first 
described by Massct & Sncllgrovc (2001). Labeling with 
subscripts 1 and 2 the inner and outer planet, respec- 
tively, in order for this mechanism to be active, the fol- 
lowing conditions must be fulfilled: 

1) the planet-to-star mass ratios (qi = Mi/M s ) must 
be such that qi > q 2 ; 

11) the separation of the semimajor axes Aa = a 2 — a\ 
must be such that Aa = 6(i?H,i + -Rh,2), where 
b < 4.5 (as we shall discuss below); 

m) q 2 must be large enough to open a gap, or partial 
gap, in the density distribution by tidal torques. 

Condition 11) above implies that Aa/a\ = (^/qi + 
v / 92)/(v^3/& — ^/(B). However, Hill stability for close 
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Figure 6. Fractional difference, Aa/ai = a2/ai — 1, between the 
semimajor axes of the interior (ai) and exterior (0,2) planets, as a 
function of the mean motion of the interior planet Hi (in scaled 
units). Different symbol sizes refer to different mean motions of 
the exterior planet, Q.2 (see the legend). The region of the graph 
below the shaded area is unstable due to planet-planet interaction. 
The Hill stability criterion adopted in this figure is given by Equa- 
tion (23) of Gladman (1993). The region above the shaded area 
does not fulfill condition 11) for gap overlap (see also Figures 7 
and 8). Resonant orbits falling in the shaded area may activate 
outward migration. 



planets on circular orbits also imposes that Aa/ai > 
2.40 + q 2 (this inequality strictly applies in the limit 
of vanishing masses and absence of gas, see Gladman 
1993), whose right-hand side is about equal to 0.26 for 
q x = M-i/Ms = 9.8 x 10~ 4 and q 2 = M 2 /M s = 2.9 x 10~ 4 , 
hence 6^2 (or 2.2, adopting a more precise determina- 
tion of the Hill stability criterion, see Gladman 1993, for 
details. See also Figure 6). Conditions i) and n) sug- 
gest that, for a Jupiter-Saturn pair, the first encoun- 
tered first-order mean motion resonance in which the 
mechanism may be activated is the 3:2 (the second-order 
5:3 commensurability being another possibility), as in- 
dicated in Figure 6. In principle, configurations exter- 
nal, but sufficiently near, to resonances may also promote 
outward migration, as we shall see in Section 3.3. It is 
worth stressing here that simple capture in a mean mo- 
tion resonance does not imply outward migration (see, 
e.g., Zhang & Zhou 2010), as we shall see in Section 7. 
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Figure 7. Surface density distribution obtained from 3D calculations of a disk with H/r = 0.04 (top) and 0.07 (middle). The angle <j>3 is 
the azimuth of the interior planet, Jupiter. The turbulence parameter at is 0.005 in both cases. A zoom of the region around the planets 
is shown in the right panels. The color scale is logarithmic and given in units of M s r^ 2 , where ri indicates the radius r = 1 (i.e., the unit 
of length). The bottom panels show the vertical stratification of the mass density, p, at the disk azimuth </>j = </>g for H/r = 0.04 (left) 
and 0.07 (right). The angle $ = w/2 — 8 is the disk's latitude and the color scale is in units of M s r~ 3 . 



Figure 7 (top and middle) illustrates the surface den- 
sity perturbed by resonant-orbit planets, derived from 
3D calculations for disks of different thicknesses (see fig- 
ure's caption). The bottom panels of the figure show the 
mass density in a vertical slice of the disk, while the plan- 
ets are aligned with the star. The exterior planet opens 
a partial gap in the case of a thinner disk, but it does so 
to a lesser extent in the other case (see bottom panels). 

The occurrence of outward migration of resonant-orbit 
planets can be intuitively understood from Figure 8, 
which reports on the results obtained from the same cal- 
culations as in Figure 7, where Saturn is caught in a 2:3 
mean motion resonance with Jupiter. The azimuthally 
averaged surface density in normalized units (long- dashed 
lines) indicates that Saturn has cleared a partial gap, 
whose inner part overlaps with the outer part of the gap 
opened by Jupiter. The short-dashed lines in the fig- 
ure represent torque density distributions (D'Angelo & 



Lubow 2008, hereafter DL08) due to Jupiter, dT /dM. 
These functions yield the total torque when integrated 
over the disk mass. The torque exerted on the interior 
planet peaks at ai ± i?n,i and is mostly comprised in a 
radial region of average width ~ 3.5i?H,i on cither side 
of the orbit (note that i?n,i > H in the case displayed in 
the left panel and R-r.i ~ H in the other case). But since 
gas depletion due to gap formation may extend somewhat 
beyond this distance (see long-dashed lines in Figure 8), 
one may allow for a maximum value of the factor b in 
condition n) above between 4 and 5. It is also important 
to notice that the orbital eccentricity of a planet acts to 
widen and smooth out gap edges (sec, e.g., Figure 2 of 
D'Angelo et al. 2006), which may also affect somewhat 
the factor b. 
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Figure 8. Azimuthally averaged surface density (long-dashed line), torque per unit disk mass (short-dashed line) exerted on the inner 
planet, and cumulative torque (Equation 20) acting on the inner planet (solid line) for the same models as in Figure 7 (cases with H/r = 0.04 
and 0.07 on the left and right, respectively). Quantity dT/dM is normalized to W 3 GM a (Mi/M s ) 2 /a\ (a\ is the semimajor axis of the 
interior planet). The surface density and cumulative torque are normalized by their absolute values at r/a\ = 2. The peaks occurring at 
the planets' orbital radii (due to mass accumulation within the Roche lobe) have been removed from the surface density profile. Results 
displayed here were obtained from 3D calculations and averaged over a few orbital periods of the outer planet. 



The solid lines in Figure 8 represent the cumulative 
torque, which is defined as 

TcmM = 2tt £^L^)r'dr'. (20) 

Looking at 7cm in the left panel of Figure 8, it appears 
clear that the positive torque exerted by the disk inte- 
rior of Jupiter's orbit is larger than that exerted by the 
disk exterior of the orbit, principally because Saturn has 
lowered the density there. The right panel illustrates the 
situation for a thicker disk, in which gas depletion op- 
erated by the exterior planet is not sufficient to reverse 
the sign of the torque. Figure 9 allows for a compar- 
ison of the cumulative torque for three different values 
of the disk thickness: H/r = 0.04 (top), 0.05 (center), 
and 0.07 (bottom) (see figure's caption for further de- 
tails). A closer inspection of dT / dM and 7cm in Figure 8 
(left) and of the cumulative torque in Figure 9 (top) in- 
dicates that the total (positive) torque is basically driven 
by Lindblad resonances and that corotation torques are 
unimportant (7cm does not vary significantly over the 
radial width of the corotation region), as also argued by 
Morbidelli & Crida (2007). 

It thus appears from Figures 8 and 9 that the discrim- 
inant factor for outward, stalled, or inward migration is 
the depth (and width) of the outer planet's gap. If the 
outer planet opened a very deep and sufficiently wide 
gap, the inner planet would be subjected only to a one- 
sided Lindblad (positive) torque exerted by the interior 
disk that, to within a factor of order unity, can be written 
as (see Lubow & Ida 2010, and referenced therein) 




Figure 9. Cumulative torque, Equation (20), exerted on the 
interior planet by a 3D disk for which the turbulence parame- 
ter is at = 0.005 and the thickness is H/r = 0.04 (top), 0.05 
(center), and 0.07 (bottom). Quantity 7cm 1s normalized to 
10 — 3 GM B 2 (Mi/M s ) 2 /ai, where a± is the semimajor axis of the 
interior planet. 



7os 



a 4 fi 2 E 



AL 



(21) 



where A = max (i7, Rn). The concept of one-sided 
torque is very useful to evaluate the presence of a tidally- 
induccd gap, i.e., condition in) above. To first-order ap- 
proximation, a density depletion begins to form when 
the one-sided torque exceeds the viscous torque (see Sec- 
tion 2.1), 7os ^ 71/, which yields a simple order-of- 
magnitude condition q 2 > 3ira t (H/a) 2 (A/a) 3 (see also 



Papaloizou & Lin 1984; Ward & Hahn 2000), or 

3/2 



Q 



>i. 



(22) 



This conditions should be regarded as a measure of how 
much the density along the planet's orbit is depleted, 
hence it can be considered a condition for gas depletion. 
A condition for tidal truncation (gap formation) is then 
j> 1 (see also Lin & Papaloizou 1986b). If predictions 
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Figure 10. Ratio of orbital frequencies of the inner to the outer 
planet showing convergence toward the 3:2 commcnsurability and 
subsequent resonant-orbit migration. The direction of migration 
of both planets is outward. The labels on the right vertical axis 
are an approximation of the semimajor axis ratio, (12/121. 

from the inequality (22) are compared with results from 
direct 3D calculations (DL08, Figures 6 and 8), one finds 
that g ~ 1 corresponds to a ~ 20% density depiction 
(relative to the unperturbed state, i.e., with no planet), 
and g « 2.7 to ~ 60% depletion. If applied to Saturn in 
the disks of Figure 8, g ss 3.4 for a density depletion of 
roughly 75% (left panels) and again g ss 1 for a density 
drop of ~ 20% (right panels). 

3.3. Orbital Migration Rates 

In order to derive migration rates, we shall assume that 
the orbits of the interior and exterior planets, Jupiter and 
Saturn respectively, are in the 3:2 mean motion resonance 
and that this resonance is maintained during migration. 
Results from a calculation that support this assumption 
are plotted in Figure 10. All calculations resulting in 
outward migration behave similarly, but there are also 
instances in which the resonance is broken (see below). 
Since the ratio 02/01 is supposed to be a constant, wc 
concentrate on the migration rate of the interior planet, 
Jupiter. We seek an expression for d of the form 

— = a rc fk 1 (E)k 2 (a)k 3 (g), (23) 
at 

where d ro f is a reference migration speed and k\ , ki , and 
&3 are dimensionless functions. To derive such expres- 
sion, we employ results from both 2D and 3D calcula- 
tions. In this section, the planets are fully formed since 
the beginning of the calculations and non-accreting. 

In Figure 11 (top-left panel), the semimajor axis evolu- 
tion of the reference model is shown for both planets (see 
figure's caption for details). The disk has an aspect ratio 
H/r = 0.04 and turbulence parameter a t = 0.005. At 
r = 1 AU, the surface density is Ei = 50gcm~ 2 at time 
t = 0. This value is chosen from the disk evolution calcu- 
lations discussed in Section 2.4, which show that > 50% 
of disks have Ei < 50 gem -2 after ~ 1 Myr. Follow- 
ing Pierens & Raymond (2011), we set 01 = 1.5 AU and 
a- 2 = 2 AU at t = (see caption of Figure 11 for further 
details), slightly outside the 3:2 mean motion resonance 
(see Figure 10). In the top-right panel of the Figure, a 
comparison between 2D and 3D calculation results is pre- 
sented. The 3D migration rate is about 25% smaller than 



the 2D one, which correction we apply to all 2D results. 
Note that this comparison is carried out at the same nu- 
merical resolution in the r <p plane (see Section 3.1). 

The parameter that may produce the largest differences 
between 2D and 3D outcomes is the disk thickness. In 
this case, however, we rely only on 3D calculations to 
approximate Equation (23). Similar plots are shown in 
the bottom panels of Figure 11, but for the orbital ec- 
centricity. For the duration of the evolution we consider, 
the eccentricity of Jupiter does not exceed ~ 0.03 in any 
of the models discussed in this section. The orbital ec- 
centricity of the exterior planet grows larger, to values 
e-i ~ 0.1 (see also Pierens & Raymond 2011). 

As mentioned above, Figures 10 and 11 indicate that 
outward migration can be activated also if orbital config- 
urations are external, but somewhat near, the 3:2 com- 
mensurability. This effect is related to the gap widths 
and the extent to which density perturbations com- 
pound. In this context, orbital eccentricity may play 
some important role, since it affects the shape of a gap 
(DAngelo et al. 2006). 

In Section 3.2, we argued that the torque exerted on 
Jupiter would tend to the one-sided Lindblad torque 7os 
(Equation 21), if Saturn (i.e., the exterior planet) carved 
a very deep and wide gap in the disk. In the opposite 
limit of very large disk thickness, H, and/or viscosity 
parameter, at, neither Jupiter nor Saturn would be ca- 
pable of depleting the disk significantly and therefore it 
is expected that the torque exercised upon the interior 
planet will be of type I (Ward 1986; Lin & Papaloizou 
1986a) 



7i 



-a 4 Q 2 E 



a \ 2 (M p 



HJ \M s 



(24) 



where, again, we neglect a factor (typically) of order 
unity in front of the right hand side. Since both 7os 
and 7i are linear in E, and since in the limit of zero 
orbital eccentricity 

dt a£lM p ' ^ ^ 

one obvious guess is to approximate k\ as a linear func- 
tion. In the left panel of Figure 12, the migration ve- 
locity a of the inner planet from calculations (symbols), 
normalized to the velocity from the reference model (Fig- 
ure 11), is plotted against the value of the normalized E. 
Here the value of the surface density is that at a distance 
of 5.5i?H,i from the (inner) planet's orbit and interior 
to it. A function proportional to E (solid line) appears 
a reasonable approximation of k\ over the range of den- 
sities shown in the Figure. Hence, we will assume that 
fci(E) = E/E re f, where the density is sampled as stated 
above. 

We note in passing that if the interior, and hence the 
exterior, planet is subjected to a torque of the type given 
in Equation (24), the resonance may be broken since the 
inner, more massive, planet may drift inward at larger 
speed than the outer, less massive, planet. This is indeed 
observed in some calculations of relatively thick disks, as 
illustrated in Figure 13 (see figure's caption for further 
details) . 

The form of function ki can be guessed following a 
similar line of argument. In their natural units of a Q 2 
(times a mass), both torques 7os (Equation 21) and 7i 



Outward Migration of Jupiter and Saturn 



13 




5000 10000 15000 20000 25000 30000 1000 2000 3000 4000 5000 6000 7000 

Time [yr] Time [yr] 

Figure 11. Left: semimajor axis (top) and eccentricity (bottom) evolution of a Jupiter-Saturn system locked in the 3:2 mean motion 
resonance, for the reference model. The planets evolve in a disk with H/r = 0.04, at = 0.005, and an initial £ at 1AU of 50gcm -2 . 
During the first 1500 years, the planets interact with each other and with the star, but do not "feel" the disk. Elliptical (osculating) orbital 
elements are obtained applying Gauss perturbation equations (e.g., Beutler 2005) in a non-rotating frame (including perturbations from 
the other planet) and then averaging over 60 and 40 orbital periods of the inner and outer planet, respectively. Right: comparison between 
evolutions of the semimajor axis (top) and eccentricity (bottom) of Jupiter, obtained for the reference model, in a 2D and 3D disk. 
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Figure 12. Normalized migration speed d/d Te { vs. the normalized (azimuthally averaged) surface density E/E re f (left) and vs. a/a Te e 
(right). The surface density is sampled at a± — 5.5-Rh.i- All the calculations are run for an evolution time between 10 and > 2 X 10 
years. To determine a, first a(t) is averaged in time over several tens of years and then a linear fit to the data is performed using a base 
time span of one to several thousand years. 



(Equation 24) scale as a 2 £ (the dependence on H/a is 
considered later), which can be seen as a measure of the 
local disk mass. Therefore, if £ was constant, a oc a? 
in units of afl (see Equation 25) and thus a would scale 
as a 3 / 2 . In the right panel of Figure 12, the migration 
speed of the inner planet from calculations (symbols), 
normalized to the reference migration speed, is plotted 



as a function of a, normalized to the semimajor axis in 
the reference model, a re f. Also in this case, the approx- 
imation seems satisfactory (solid line) and so we shall 
assume that k^a) = (a/a ro f) 3 / 2 . 

In Section 3.2, it was anticipated that the depth and 
width of the density depletion produced by the exterior 
planet play a fundamental role in determining magni- 
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Time [yr] 

Figure 13. Scmimajor axes vs. time of a Jupiter-Saturn type 
system evolving in a 3D disk with H/r = 0.07, at = 0.005, and 
an initial £ at 1 AU of 50 gem -2 . Each semimajor axis is normal- 
ized to its initial value. The planets are subjected to disk torques 
after the first 1500 years. The planets undergo divergent inward 
migration. 
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Figure 14. Normalized migration speed a/a re f for various values 
of the ratio <?/g r ef, calculated for the exterior planet, according to 
Equation (27). All models are run for an evolution time of between 
6 x 10 3 and 1.2 x 10 4 years. The thick horizontal line segment 
represents the type I migration (see Equation 24) that the inner 
planet would be subjected to if the disk had a relative thickness 
H/r = 0.1. 

tudc and direction of the interior planet's migration, and 
hence of the pair as a whole. Quantity g (Equation 22) 
can be used as a proxy to discriminate among the vari- 
ous situations, i.e., different combinations of H/r and at 
(and q) that may affect E. On account of the compact 
orbital configuration, since q\ > qi we have that g\ > gi 
(unless H/r and/or at are rapidly varying with disk ra- 
dius). If gi ^> 1, then g\ 3> 1, and the interior planet is 
likely subjected to a torque whose limit is the one-sided 
Lindblad torque in Equation (21) and the two planets 
migrate outward. Otherwise, if gi <C 1, then g2 -C 1, 
and migration will be dictated by a type I torque (Equa- 
tion 24) and be directed inward. 

The migration velocity of the inner planet should de- 
pend on both g\ and g2- However, in the present context 
<?2 should have the larger impact of the two and there- 
fore, for the sake of simplicity, we assume that function 
k% depends only on g2- In particular, referring to the 



value of 92 in the reference model (Figure 11) as g re { and 
indicating 52 simply as j, £3 will be approximated as 
a function of g/g Te i- Since it is unclear how the agree- 
ment between 2D and 3D calculations varies as a func- 
tion of the disk thickness (it should likely worsen as H 
increases), we only use 3D calculations to find an approx- 
imation to function k^. Figure 14 shows the ratio d/d ro f 
obtained from calculations for various values of the ra- 
tio g/ g-cd- The thick horizontal line in the plot indicates 
the type I migration speed that would apply to the inner 
planet if H/r = 0.1, corresponding to the value used for 
left-most data point on the graph. For reference, the nor- 
malized migration speed corresponding to the one-sided 
torque in Equation (21) would be ~ 90. The broken line 
is a linear interpolation of the data, which will be used 
as a representation of fca . 

3.4. Approximation of the 3:2 Resonant- orbit 
Migration Velocity 

Summarizing the results of Section 3.3, we write the 
migration speed of the interior planet as 

do f E \ f a \ 

^ =dr H^Jw Wfef) ' (26) 

where the dimensionless function k^ is obtained via lin- 
ear interpolation of the numerical data in Figure 14, i.e., 
the thick solid line in the Figure. The reference values 
in Equation (26) are taken from the reference model, 
corrected for 3D effects, as discussed in Section 3.3: 
d rcf = 2.7 x 10" 6 AUyr- 1 for E ref = 42gcm~ 2 and 
a ro f = 1.57 AU. Recall that both E and E re f are evalu- 
ated at a distance of 5.5 i?H,i interior of the inner planet's 
orbit. The argument of function £3 is given by 

9 _ /at.ref ( H le { \ f #H \ 

-g-f-y—y^r) [tj ' { 7) 

where a tj ref = 0.005-^/1 AU/a rcf , (H/a) Ic{ = 0.04, and 
A = max (H, Rn). Recall that here parameters g and g rc { 
are computed from Equation (22) applied to the exterior 
planet. For a mass ratio q different from that of the refer- 
ence model, the right-hand side of Equation (27) should 
be multiplied by q/q Ye f. Equation (26), without the cor- 
rection due to 3D effects, is also in reasonable agreement 
with the results presented by Pierens & Raymond (2011) 
in their Figure 21, for a disk of 0.4 Mj within 1.5 AU 
(E ~ 500gcm~ 2 at 1 AU). As explained above, the exte- 
rior planet's orbit may not always be resonant with that 
of the interior planet, when migration is inward (see Fig- 
ure 13). Nonetheless, for the outer planet's orbit we set 
a 2 = (3/2) 3 /3 ai . 

Equation (26) predicts stalling points (where &3, and 
hence d, is ps 0) at <?o ~ 0.8g le {, which will be regarded as 
a nominal value. But notice that since the approximation 
to &3 is sampled at a limited number of points, its zero 
could be located at a somewhat different abscissa, yet 
it is located between 0.7c/ re f and g re {. In principle, there 
could be multiple stalling points in a disk (if disk proper- 
ties are not monotonic), which would represent locations 
of stable equilibrium since they are convergent radii for 
the planets' semimajor axes. The argument, however, is 
based on the assumption that the 3:2 commcnsurability 
is preserved even for inward migration. This is not true 
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Figure 15. Ratio g/g Te { (Equation 27) as a function of the ex- 
terior (bottom axis) and interior (top axis) planet's orbital radius. 
Here, H/r is approximated by Equation (16). The shaded area 
indicates the region where di changes sign: di > (di < 0) above 
(below) the shaded area. The horizontal line within the boundaries 
of the shaded area is the nominal value for stalling migration (see 
the text) . The various curves represent g/g IC { for different "slope" 
parameters, j3 (see the legend). Thin and thick curves correspond 
to v\ = 4 X 10 -6 and 1.6 X 10 -5 r\ Oj, respectively. 

in general (as illustrated, for example, in Figure 13), but 
it further requires that the inward migration of the ex- 
terior planet be faster than that of the interior planet 
and that the condition for gap overlap (see condition u) 
in Section 3.2) be satisfied. When these conditions are 
not met, the pair of planets can migrate past a stalling 
point, toward the star. This may happens, for example, 
if disk temperature and viscosity are low enough so that 
both planets open up a gap and drift according to type II 
migration before capture into resonance occurs. 

Some constraints on the range of outward migration 
predicted by Equation (26) can be derived for the disk 
models discussed in Section 2.4. The disk thickness, H, is 
affected by internal (viscous) heating typically over the 
first few million years of evolution (see Figure 3, left). 
If we neglect that source of heating, H can be approx- 
imated by using Equation (16) and then Equation (27) 
can provide a rough measure of a disk's radial range over 
which the planets can drift away from the star, as shown 
in Figure 15. The addition of internal heating would raise 
the value of H, hence reducing the ratio g/<? ro f, which 
suggests that outward migration may not be activated 
in the warm interiors of a young disk. The radial region 
over which a curve extends above or, to some degree, in- 
side the shaded area is favorable to outward migration. 
The intersection of a curve with the horizontal line in the 
shaded area gives the nominal radius of the stalling point 
of each planet (see above). According to the plot, out- 
ward migration of a Jupiter- mass planet locked in the 3:2 
mean motion resonance with a Saturn-mass planet can- 
not proceed beyond ~ 7 AU (top x-axis). The numerical 
experiments discussed in Section 4 are broadly consistent 
with this prediction. 

The validity of Equation (26) for planet-to-star mass 
ratios different from those adopted here (q\ = Mi/M s = 
9.8 x 1CT 4 and q 2 = M 2 /M s = 2.9 x 10~ 4 ) was not inves- 
tigated. Parameter g (Equation 22) is oc q 2 if i?H < H 
and oc y/q~2 otherwise. Therefore, a faster outward mi- 
gration might be expected as q 2 increases, provided that 
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Figure 16. Shaded area represents orbital frequency ratios, as a 
function of the interior planet mass, which may possibly activate 
outward migration of a pair of planets whose mass ratio is such that 
9i/<72 = M1/M2 = 3. The upper boundary of the shaded area is 
given by condition 11) of Section 3.2. The lower boundary is given 
by a Hill stability criterion for circular orbit planets (Equation 23 
of Gladman 1993), in absence of gas. It is assumed that parameter 
g, defined in Equation (22) and applied to the exterior planet, is 
large enough to allow for sufficient gas depletion. The area bounded 
by the thicker lines represents this condition in terms of the ratio 
aila\ (labeled on the right vertical axis). 

the ratio qi/q 2 stays roughly constant. Condition n) of 
Section 3.2, along with the requirement of Hill stability, 
suggests that, as q 2 increases beyond ~ 0.001, the 2:1, 
rather than the 3:2, commcnsurability may be available 
to activate outward migration for a ratio q\/q 2 ~ 3, in 
accord with the findings of Crida et al. (2009). This is 
schematically illustrated in Figure 16 (see figure's cap- 
tion for further details). As the ratio q\/q 2 approaches 
1, the shaded area in the graph shifts upward. 

4. LONG-TERM 3:2 RESONANT-ORBIT 
MIGRATION OF JUPITER AND SATURN 

The results from the disk models of Section 2.4 can 
be combined with the results of Section 3.4 to study the 
migration of a resonant-orbit pair of planets with masses 
corresponding to those of Jupiter and Saturn. Here we 
shall assume that, by a time r p , both planets have fully 
formed, i.e., they have reached their final masses, and 
their orbits have become locked in the 3:2 mean mo- 
tion resonance. Planet formation calculations of a giant 
planet via core nucleated accretion (e.g., Hubickyj et al. 
2005; Alibert et al. 2005; Lissauer et al. 2009; Movshovitz 
et al. 2010; Mordasini et al. 2011) indicate that the for- 
mation time of Jupiter is greater than ~ 1 Myr, although 
this timescale is affected by the formation of the solid 
core and thus by the distance where the core forms. The 
formation of Saturn, the exterior planet, seemingly takes 
somewhat longer (see, e.g., Pollack et al. 1996; Dodson- 
Robinson et al. 2008; Benvenuto et al. 2009). In addition 
to the formation time, there is the time required by con- 
vergent orbital migration to bring the planets into mean 
motion resonance (which is also included in r p ). Assum- 
ing that t p is determined mainly by the formation time 
and for lack of better constraints, we choose three refer- 
ence times t p of 1, 2 and 3 Myr (e.g., Kenyon & Bromley 
2009; Bromley & Kenyon 2011). 

It is important to bear in mind that longer timescales 
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are possible, whereas it is unclear if shorter timescalcs arc 
feasible. In fact, there are also observational constraints 
suggesting that Jupiter formed after several million years 
(see Scott 2006, and references therein). 

The disk evolution models presented in Section 2.4 
are recalculated, starting from time t p and using Equa- 
tions (26) and (27), to integrate the orbital radius of 
the interior planet. The orbit of the exterior planet is 
constrained by the 2:3 resonant-orbit requirement with 
the interior planet (see Section 3.4). In particular, the 
disk models provide the quantities E, at, and H , used in 
Equations (26) and (27), as they evolve over time. The 
effects of the torques exerted by the planets on the disk 
are not taken into account in the ID models because, for 
current purposes, the feedback of the tidal field on the 
disk and its effects on the migration rates are included in 
the 2D and 3D calculations discussed above. The mass 
of both planets is constant, but we will contemplate the 
impact of gas accretion in Section 5. Since t p > lMyr, 
useful disk models are those for which td > 1 Myr (and 
< 20 Myr). 

The evolution of a\, the interior planet's orbital radius, 
is illustrated in Figure 17 for selected cases (see figure's 
caption for details). The complete list of the asymptotic 
values of a±, Oqc, is reported in Table 2 for r p = 1 Myr. In 
some cases, the orbital radius remains nearly unchanged. 
This happens as a result of E (around r ~ a\) being too 
low, at time t = 1 Myr, for any significant amount of 
angular momentum to be transferred to/from the disk. 
Obviously, this result holds for t p > 1 Myr. 

The outcomes of Figure 17 and Table 2 can be inter- 
preted with the aid of Figure 15. For the highest viscos- 
ity regime, outward migration is not activated, and both 
planets migrate inward regardless of other disk parame- 
ters. At the intermediate viscosity, cti may be positive, 
but the nominal stalling radius is within 2 AU, hence the 
interior planets may not proceed beyond this distance. 
The lowest viscosity regime offers the widest range of 
outward migration, resulting in nominal stalling radii of 
« 6.5AU = 1/2), « 2.4AU {fi = 1), and « 1.6AU 
(fi = 3/2). Because of the variation of v with radius, as 
the density steepens, the nominal stalling radius moves 
inward. At even smaller viscosity, migration of the res- 
onant pair may proceed to larger distances, but in this 
case there are at least two possible issues that may arise. 
One is related to the disk lifetime, which increases as 
v decreases (see Table 1). The other is related to the 
mode of migration of the exterior planet prior to reso- 
nance capture, which may transition to type II at low 
enough viscosity, hence convergent migration toward the 
interior planet may be compromised (see Section 7). 

Figure 15 can also assist in extending the results il- 
lustrated in Figure 17, and reported Table 2, to differ- 
ent initial orbital radii. Assuming that viscous heat- 
ing does not represent a major source term in the en- 
ergy budget, Equation (14), a pair of planets that be- 
come locked in the 3:2 mean motion resonance inside the 
stalling radii (of each planet) will migrate toward those 
locations. Whether or not the planets may reach those 
radii depends on the gas density level in the disk. At 
v\ = 4 x 10~ 6 rffii, models for which aoo sa 6.5 AU 
for P = 1/2, or « 2.4 AU for /3 = 1, or w 1.6 AU for 
[3 = 3/2, have reached their stalling radii. If the pair 
becomes locked into resonance outside the stalling radii, 



the planets will converge toward, or transit across, them 
depending on the disk conditions (sec discussion in Sec- 
tion 3.4). 

Among the sets of parameters listed in Table 2, only 
five are compatible with the outward migration of the 
interior planet beyond ~ 5 AU (when resonance lock- 
ing occurs at ~ 1 AU) . Only two sets of parameters 
remain compatible with this requirement if the forma- 
tion timescale is r p = 2 Myr, and only one if t p = 3 Myr, 
as illustrated in the left panel of Figure 18 (see figure's 
caption for further details). 

In order to derive a distribution of the asymptotic or- 
bital radii of the interior planet, a^, for each pair of 
values (Mp,/3) listed in Table 2, v\ and /41 are var- 
ied randomly between the corresponding minimum and 
maximum values reported in the table, and r p is varied 
randomly between 1 and 3 Myr. A total of over 1200 
models were computed and the histogram of the results 
is shown in the right panel of Figure 18. Overall, there 
is a 97% probability that the interior planet will achieve 
an asymptotic radius doo < 3 AU and a 98% probability 
that doo < 4 AU. 

5. GAS ACCRETION AND PLANET GROWTH 

As mentioned in Section 3.2, the first condition nec- 
essary to activate the outward migration mechanism of 
resonant-orbit planets is that the interior planet's mass 
must exceed that of the exterior planet. In the limit 
of equal mass planets, one expects the outer Lindblad 
torque exerted on the exterior planet to overcome the in- 
ner Lindblad torque exerted on the interior planet (see, 
e.g., Morbidelli & Crida 2007). Several outcomes are 
then possible, including breaking of the resonance, scat- 
tering, and inward type II migration of both planets. 

The neglect of gas accretion, especially on the exterior 
planet, represents possibly the most serious limitation of 
this mechanism. Hydrodynamical calculations can pro- 
vide maximum, or disk-limited, gas accretion rates for 
such planets. Although they do not necessarily represent 
the actual accretion rates, formation models of Jupiter 
(Lissauer et al. 2009) indicate that, once runaway ac- 
cretion begins, a giant planet does grow at a disk-limited 
accretion rate. In a disk with H / r < 0.05 and at < 0.005, 
an isolated Saturn-mass planet may accrete gas at a rate 
a few times as large as that of a Jupiter-mass planet at 
the same location in the disk (see Lissauer et al. 2009, 
and references therein). In a disk with H/r 3> 0.05 or 
at ^s> 0.005, these rates would be comparable (DL08). 
Even assuming the same accretion rate for both plan- 
ets, M p , the initial growth time, M p /M p , of Saturn is 
shorter than that of Jupiter hence Saturn may approach 
the mass of Jupiter more or less quickly, depending on 
the local values of E and H/r. Furthermore, there is no 
obvious reason as to why Jupiter and Saturn should ac- 
crete gas at a disk-limited rate and then suddenly stop 
accreting despite the continuing supply of gas from the 
disk. 

The evolution of one high-resolution 3D calculation 
was continued by allowing the two planets to accrete gas 
following the procedure outlined in DL08, modified to 
account for the different local dynamical times 6 (i.e., the 
timescale for mass removal depends on the planet's or- 

6 If accretion proceeds through a disk around the planet, as 
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Figure 17. Evolution of the orbital radius (a±) of the interior planet computed by integrating Equation (26) along with the disk models 
of Section 2.4 and assuming a formation timescale t p = lMyr (see the text). The initial disk mass, in units of M s , is indicated in the 
top-right corner of the left panels. The initial "slope" parameter of S is /3 = 1/2 (left), 1 (center), and 3/2 (right). The curves in each 
panel refer to different values of v\ and /41 (see Table 2 to identify each curve). A complete list of the asymptotic values, <2oo, is reported 
in Table 2. It is assumed that a 2 = (3/2) 2 / 3 ai. 
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Figure 18. Left: orbital radius of the interior planet vs. time, assuming a formation timescale r p = 2 Myr and 3 Myr. Thin and thick 
curves represent, respectively, models with an initial disk mass MS = 0.044 (/41 = 10) and 0.088 M s (/41 = 100). In all cases, parameters 
are /3 = 1/2 and v\ = 4 X 10 _6 r^f2i. For r p = 2 and 3 Myr, these are the models that show the longest range of outward migration 
among those listed in Table 2. Right: histogram of the asymptotic orbital radius of the interior planet obtained by varying randomly the 
parameters v\, /41, and r v . See the text for further details. 



suggested by the bottom-left panel of Figure 7, M p is equal to the 



rate at which the nebula supplies this disk, hence the connection 
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Table 2 

Asymptotic Orbital Radii for t p = 1 Myr 
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a Asymptotic value of the semimajor axis, in AU, of the interior planet's orbit. 
b Initial "slope" of the disk's surface density. 

c Kinematic viscosity at ri = 1 AU in units of r\ f2i = (G M s ri) 1 / 2 . 

d Rate of EUV ionizing photons emitted by the star in units of 10 41 s — 1 (see Equation 5). 
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Figure 19. Disk-limited gas accretion rates, in units of M s yr -1 , 
of the interior (thin curve) and exterior (thick curve) planets. The 
pair is locked in the 3:2 mean motion resonance and the year count 
starts after about 6700 years of evolution. The plot shows only a 
small time interval to highlight the accretion modulation. The nu- 
merical resolution is such that there are ~ 27 3 and ~ 20 3 grid cells 
in the Hill sphere of the interior and exterior planets, respectively. 

bital frequency). Disk conditions are similar to those 
applied to the reference model of Section 3.3. The disk- 
limited accretion rates of the two planets are shown in 
Figure 19, where the thicker curve refers to the exte- 
rior planet. Both accretion rates are modulated (notice 
the logarithmic scale) over the orbital period due to the 
eccentric orbits (D'Angelo et al. 2006), although addi- 
tional modulations may be present due to the resonant 
forcing. For conditions simulated in Figure 19, the in- 
tegrated values of M v are very similar, differing by less 
than 10%. These rates would yield a growth time for 
the exterior planet of a few times 10 4 years, typically 
much shorter than the migration time (see Figure 17). 
For comparison, the growth time of the interior planet 
would be ~ 10 5 years. The net effect of gas accretion 
could produce a mass ratio M2/M1 close to 1, or possi- 
bly larger (since gas starvation would likely occur for the 
interior planet first, see Figure 3). 

with the orbital frequency. 



The orbital evolution of the accreting planets, moni- 
tored over ~ 1000 years, does not show any significant 
deviation from the evolution of the non-accreting planets 
(the planet masses are fixed in this case, on account of the 
little variations expected over that timescale). However, 
as explained in Section 6, one or more accreting planets 
may change the steady-state structure of the inner disk, 
affecting the migration behavior of the planets. 

6. GAS ACCRETION AND EFFECTS ON THE DISK 

Resolving the problem of the rapid growth of the ex- 
terior planet, in a still relatively massive disk, would re- 
move only one issue posed by gas accretion. In fact, even 
if Saturn suddenly stopped accreting, gas accretion onto 
Jupiter would continue to pose a problem. The issue here 
is not related to the growth of the interior planet's mass, 
but rather to the modification of the mass flux through 
the disk, across the planet's orbit. 

As explained in Section 2.1, in a stationary disk the ac- 
cretion rate is M = 3ttvH and nearly independent of the 
radius r. Accretion on the interior (or exterior) planet 
would change M. This phenomenon was analyzed in de- 
tail by Lubow & D'Angelo 2006 (hereafter LD06) for the 
case of a single planet. The generalization to two planets 
can be performed by introducing an average accretion ef- 
ficiency that quantifies the amount of gas accretion onto 
both planets, relative to an average local accretion rate 
through the disk . However, here we wish to consider 
the situation in which the interior planet accretes gas, 
but the exterior planet does not. Therefore, the formal- 
ism of LD06 can be applied in a straightforward manner. 
Let us indicate with M c = M = "invY, the accretion rate 
sufficiently far from the exterior planet's orbit (so that it 
is basically unperturbed) and with M; the accretion rate 
inside the orbit of the interior planet. If there was no 
sink in the disk, then M; = M e . Yet, since some mate- 
rial is removed by the planet, in general Mj < M e , which 
implies a reduction of £ in the inner disk with respect 
to the same disk without the planet. This reduction de- 
pends on the planet's accretion efficiency £, defined as 

7 As pointed out by LD06, for a given disk, there is a planet mass 
(Mi + M2, in this case) beyond which the couple exerted by the 
planct(s) will make the accretion disk evolve toward a decretion 
disk (Pringle 1991). 
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the ratio of M p to the accretion rate interpolated at the 
planet's orbital radius. Since Mi = M c — M p , one finds 
that Mi = M / (£ + 1). This result formally applies if the 



disk's inner boundary, 



is much smaller than a\. In 



general, one finds that (see LD06) 



M = 



1 + (1 - \Amin/ai)£ 



(28) 



According to the Equation (28), the surface density in 
the inner disk is then expected to be reduced by a fac- 
tor of order £ + 1 (assuming that v remains unchanged), 
relative to the situation in which M p = 0. Accordingly, 
the (positive) Lindblad torque foe E, see Equation 21) 
exerted by the inner disk on the (inner) planet is also ex- 
pected to decrease. Such effect may slow down outward 
migration and may even tip the torque balance in favor 
of the negative torque acting on the interior planet. 

We estimate the efficiency of accretion, £, allowing 
for accretion on the interior planet only by applying the 
steady-state solution given by Equation (19) of LD06 as 
initial condition and using the iteration procedure out- 
lined in LD06. The disk configuration is as that of the 
reference model in Section 3.3, except for a somewhat 
smaller disk's inner radius of 0.15 ai and for the initial 
location of the exterior planet (02/01 = 1.315, see Fi- 
grue 10). We find that £ ~ 6. There are fluctuations over 
time of both M p and M (see Figure 19) and therefore the 
accretion efficiency is taken as the ratio of averaged quan- 
tities. The estimated variation is A£ ~ 1. The accretion 
rate ratio is Mi/M c ~ 0.2, whereas the corrected value, 
that is the ratio estimated for r min /ai < 1, is ~ 0.14 
(see Equation 28). 

In Figure 20 (left panel), the surface density after 1000 
orbits of the interior planet (solid line) is compared with 
the steady-state solution (Equation 19) of LD06 for a disk 
with a single planet (dashed line). The reduced mass ac- 
cretion past the planets, Mi, results in a lower surface 
density (compare with the long-dashed line in the left 
panel of Figure 8). The migration of the resonant-orbit 
pair (initially placed in the 3:2 mean motion resonance) 
in the stationary surface density of Figure 20 (left panel) 
is shown in the right panel. The reduced positive Lind- 
blad torque by the inner disk cannot overcome the nega- 
tive torque by the disk outside the planet's orbit, result- 
ing in a migration speed &i < (note that the units of 
time in Figure 20 are initial orbits of the inner planet, 
not years). 

7. THE 2:1 MEAN MOTION RESONANCE 

A pair of planets undergoing convergent migration will 
first cross the 2:1 commensur ability, before approach- 
ing the 3:2 mean motion resonance. While the latter 
resonant-orbit configuration may activate the outward 
migration mechanism discussed here (if M1/M2 ~ 3), 
the former may not (see Figures 6 and 16). But a pair 
of giant planets interacting with a gaseous disk can in- 
deed be caught in this resonance, as shown by several 
studies (see, e.g., Kley 2003; Kley et al. 2004; Pierens & 
Nelson 2008; Zhang & Zhou 2010). We therefore seek 
conditions such that the outer planet may, or may not, 
overcome the barrier represented by capture in the 1:2 
mean motion resonance, while migrating toward the in- 
ner planet. 



In the case of convergent migration, the condition for 
capture of the exterior planet in a resonant orbit with 
the interior planet requires that the relative migration 
speed be such that 



da IC 



dt 



< 



T, 



(29) 



where d rc i = a<z — <ii, Aa rGS is the resonance amplitude in 
scmimajor axis, and 7} is the resonant libration period, 
i.e., the period of the critical angle (e.g., Michtchenko 
et al. 2008) Vi = 2(.M 2 +tu 2 ) - (M\+mi) -mi, where 
Ai indicates the mean anomaly and w the argument of 
periapsis. Recall that, in the case of a more massive inte- 
rior planet, ipi is a real resonant angle, as it displays dy- 
namical libration. For low-eccentricity orbits (e < 0.1), 
Mustill & Wyatt (2011) (see also Quillen 2006) approxi- 
mate the critical relative velocity for capture of the outer 
planet in the 1:2 mean motion resonance as 



da rc \ 


< 1.2 


'MA 


dt 







4/3 



(30) 



where it is assumed that f2x/^2 = 2. Capture appears 
to be probabilistic at higher eccentricities, with a critical 
relative velocity that becomes somewhat larger for larger 
eccentricities. Here we should stress, however, that the 
estimate of the critical velocity in Equation (30) assumes 
captures of "particles", i.e., Mi 3> M 2 . If the planet 
transits the 1:2 orbital resonance with the interior planet, 
capture in the next first-order resonance, the 2:3, requires 
that I Orel I < U(M 1 /M s ) i / 3 a 2 n 2 . 

We shall assume that the inward migration speed of 
the interior planet is negligible compared to that of the 
exterior planet, thus ct ro i Rj a 2 . The exterior planet 
may in principle undergo a mode of rapid migration 
dominated by corotation torques (type III). D'Angelo 
& Lubow (2008) separated this mode of migration from 
type I mode by analyzing the torque density distributions 
and the fluid trajectories in the corotation region of the 
planet. They found that two conditions must be satis- 
fied for the activation of type III migration. The first is 
that the migration timescale across the coorbital region 
is shorter than the timescale required to clear a gap over 
that same region. The second is that the unperturbed 
surface density at the planet location is such that 



q 2 £ 

m7 



> 



(31) 



where disk quantities are sampled at a = a 2 8 . If the first 
condition is fulfilled, the second condition requires that 
E > 10~ 3 Ms a~ 2 for H/a > 0.03, which corresponds 
to a density in excess of 10 gem" 2 in the disk region 
within 3AU of the star. According to Figures 4 and 5, 
this requirement is not met, thus we can assume that 
the outer planet migrates at a rate in between type I 
and type II migration rates. Assuming a speed of order 

8 The condition represented by Equation (31) is also consistent 
with the case reported by Masset & Snellgrove (2001), which show 
an exterior planet migration dominated by corotation torques. In 
that case, a\s/M s = 6 x 10~ 4 , H/r = 0.04, and a 2 = 2 (in their 
units). Hence, evaluating at the initial position of the exterior 
planet, one has a|S/M s = 2.4 X 10~ 3 > {H/a 2 ) 2 . 
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Figure 20. Left: azimuthally averaged surface density (solid line), in units of M B a^ 2 , of a disk with a pair of planets, initially placed 
in the 3:2 commcnsurability. The exterior planet is non-accreting. The interior planet accretes gas with an efficiency parameter £ ~ 6 
(see the text). The dashed line indicates the steady-state solution of LD06 (with a single planet) for £ = 6 and an inner disk radius of 
0.15 ai. Right: migration track of the interior planet. Data are averaged over ~ 50 orbital periods. The planet begins migrating after 1000 
(initial) orbits. The depiction of the inner disk (compare with Figure 8, long-dashed line) is sufficient to deactivate the outward migration 
mechanism (compare with Figure 11. top, which has different units on the time axis). 



type I and using Equation (24), one finds 



da 



rcl 



dt 



(32) 



all disk quantities being evaluated at a = a 2 - In the equa- 
tion above, a numerical factor of order unity multiplying 
the right-hand side is neglected. This is done to account 
for the fact that the density is partly depleted and the 
actual migration rate deviates somewhat from the type I 
rate (see DL08). Moreover, we find that Equation (32) 
gives a reasonable order-of-magnitude approximation to 
numerical results. 

Therefore, capture of the exterior planet in the 1:2 
mean motion resonance with the interior planet may oc- 
cur if the unperturbed surface density is lower than a 
critical value, so that 
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If H/a > 0.03, then a surface density £ < 3 x 
10 -4 A/ S a~ 2 , or about 650gcm~ 2 at m 2 AU, may be 
sufficient to allow for capture in this resonance. If the 
exterior planet crossed the 1:2 commcnsurability while 
it had a much lower mass, say ~ 20 Earth masses, this 
critical density would not decrease, even accounting for 
a numerical factor of 4-5 in Equation (32) and restoring 
a full type I migration. 

The inequality in Equation (33) suggests that if 
a|£/M s > 5 x 10" 4 in a disk with H/r ~ 0.04, the ex- 
terior planet may transit the 1:2 mean motion resonance 
with the interior planet. These conditions are realized, 
for example, in the model of Masset & Snellgrove (2001) 
(see their Figure 1) and in the models of Pierens & Nelson 
(2008) (see their Figure 5), considering that the density 
must be rescaled at the radius of the resonance crossing, 
where a 2 ~ 1.5 (in the units of Masset & Snellgrove 2001) 
and a 2 ~ 1.3 (in the units of Pierens & Nelson 2008). 
The solar nebula models considered here, however, sug- 
gest that if the orbits approach the 1:2 commcnsurability 
in the inner disk, after ~ 1 Myr, then capture is likely in 
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Figure 21. Ratio of the mean motions vs. time for a pair of 
planets undergoing convergent migration. The value of the disk's 
surface density at 1 AU is indicated at the bottom of each panel, in 
units of g cm -2 , for the curves of different thickness. The top panel 
illustrates evolutions for typical densities obtained from the disk 
models of Section 2.4. Circles are predictions from Equation (34). 
The inset shows the migration tracks of the interior planet after 
resonance locking occurs. At higher densities, capture of the exte- 
rior planet transitions from the 1:2 to the 2:3 commcnsurability, as 
illustrated in the bottom panel. 



a statistical sense (see Figure 5). 
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Equation (33) represents only an approximate condi- 
tion for resonance locking, because orbital eccentricity 
may play some role and because Equation (30) was not 
derived for the capture of similar mass bodies. In or- 
der to provide a further test on conditions that may lead 
to locking of the exterior planet in the 1:2 commensura- 
bility, calculations along the lines of those presented in 
Section 3.3 are performed for varying initial surface den- 
sity at 1 AU, Ei. The disk conditions are those used for 
the reference model in Figure 11, except that the exterior 
planet is placed initially on a circular orbit at a distance 
of 2.8 AU from the star, outside the resonance location 
with the interior planet, and both begin migrating after 
500 years. Results from these calculations are illustrated 
in Figure 21, which shows the ratio of the mean mo- 
tions. The top panel refers to values of Ei compatible 
with those found in the disk evolution models at the time 
of planet formation (see Figure 5, left). In these cases, 
the exterior planet becomes locked in the 1:2 orbital res- 
onance. Since £1 = — (3/2)Od/o, assuming that £li is 
nearly constant, then 

d /O x \ 3 fflA ( a \2 f M 2 \ /a 2 E\ , N 

*U)~-5(n;)(if) UXmtH (m » 

where quantities depending on a are evaluated at a = 
CL2- Predictions from Equation (34) are superimposed 
(circles) to f2i/f22 curves in Figure 21, indicating that 
the exterior planet does approach the interior planet, at 
least initially, with a radial speed on the order of the 
relative velocity given by Equation (32). 

The bottom panel of Figure 21 shows cases with higher 
initial densities in which capture of the exterior planet is 
in the 1:2 or the 2:3 orbital resonance. There is overall 
agreement with Equation (33), and transit across the 
1:2 orbital resonance is obtained for r 2 Ei/Af s = 6 x 
10 -4 . We did not investigate the exact density value at 
which locking transitions from one to the other resonant 
configuration, but it is likely that there exists an interval 
of values for which the result is stochastic. 

The inset in the top panel of Figure 21 illustrates the 
migration of the interior planet after resonance locking. 
In the bottom panel, migration is outward for the case 
that shows locking of the exterior planet in the 2:3 mean 
motion resonance with the interior planet, inward in the 
other cases. These results confirm what was argued 
above and suggested by Figures 6 and 16: the 2:3 orbital 
resonance leads to outward migration, the 1:2 resonance 
does not. 

A small disk aspect ratio may help preventing capture 
of the exterior planet in the 1:2 commensurability with 
the interior planet. However, the condition for gap for- 
mation (Equation 22) suggests that the migration of the 
exterior planet may transition to type II, and hi cannot 
be neglected. In this case, since v cx a 13 and a ~ vja oc 
a* 3-1 , convergent migration requires that j3 > 1. By 
applying Equation (30), one finds that if the kinematic 
viscosity at 1AU is v x < (A/i/Af s ) 4 / 3 (l AU/a 2 ) (/3 - 1/2) 
then convergent type II migration may still lead to or- 
bital locking in the 1:2 mean motion resonance. However, 
this estimate is complicated by the fact that if the planet 
mass becomes larger than the local disk mass, a likely sit- 
uation at late evolutionary times, a is also proportional 
to a 2 T,/M p due to intervening inertia effects (see, e.g., 



Syer & Clarke 1995; Ivanov et al. 1999). Hence, the in- 
terior planet would likely slow down before the exterior 
planet would. 

7.1. Saturn Formation within the 1:2 Orbital 
Resonance with Jupiter 

We shall consider here the possibility that Saturn forms 
within the 1:2, but outside the 2:3, commensurability 
with Jupiter, while both planets are beyond several AU 
from the Sun. If during the course of its evolution Sat- 
urn remains inside the 1:2 mean motion resonance with 
Jupiter, then capture in the 2:3 orbital resonance is still 
possible. 

In order to maintain such a compact orbital configu- 
ration throughout the evolution of the two planets, the 
migration rates must be very similar over time. In fact, if 
the constraint 1.31 < 02/01 < 1.59 has to be preserved, 
a change of this ratio of at most 20% over the forma- 
tion timescales basically implies that a 2 ja\ is roughly 
constant. Therefore, by taking the time derivative, one 



But since the interior planet has to grow faster than the 
exterior planet does, their migration rates are bound to 
differ, at some point in time at least. Thus, the condition 
in Equation (35) is unlikely to be (always) satisfied and 
the difference a 2 — a\ will either increase or decrease. 

If Saturn forms within the 1:2 orbital resonance and 
Jupiter has still to acquire most of its mass, there will be 
a phase when the migration rate of Jupiter significantly 
exceeds that of Saturn (presumably, around the time of 
runaway gas accretion; see DL08). Hence, it is most 
likely that Saturn is left behind and becomes probably 
trapped in the 1:2 commensurability. Instead, if Jupiter 
has already acquired the bulk of its mass, the opposite is 
likely to occur and most probably Saturn becomes locked 
in the 2:3 orbital resonance while it is still growing. 

We consider this last situation in some details by per- 
forming 3D calculations similar to those of the reference 
model of Section 3.3, but applying a lower mass to the 
exterior planet: q = M 2 /M s — 10 -4 and 2 x 10~ 4 . 
Equation (27), appropriately modified for mass ratios 
Q 7^ Qicf (see discussion after Equation 27), suggests 
that if q ~ 10~ 4 migration of the pair is inward, but 
if q w 2 x 10 -4 then migration is outward. Direct calcu- 
lations agree with these predictions, as indicated by the 
cumulative torques shown in the top-left panel of Fig- 
ure 22 for three different values of M 2 (sec caption for 
details) . The density depiction due to the tidal torques of 
the exterior planet amounts to ~ 35% for M 2 = 10 -4 M s 
and to ~ 60% for M 2 = 2 x 10 _4 M S (top-right panel), 
in accord with the expectations of Equation (22). The 
bottom panels of Figure 22 show vertical distributions of 
the mass density (see figure's caption). 

Therefore, if Saturn grows while locked in the 2:3 or- 
bital resonance with Jupiter, there is a mass smaller 
than Saturn's final mass for which migration stalls and 
then reverses, as suggested by 7cm as function of M 2 
in Figure 22 (top-left panel). This situation would pre- 
vent Jupiter from reaching the inner disk regions, unless 
Saturn achieved the mass for migration reversal when 
Jupiter is already there. For H/r « 0.04, such mass 
is between 10 _4 M S and 2 x 10~ 4 M S (and presumably 
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Figure 22. Top: cumulative torques (left) exerted on the interior planet and averaged surface density (right), obtained from 3D calculations 
whose parameters are as in the reference model of Section 3.3. Different lines types indicate different masses of the exterior planet (M2/ M s ): 
2.9 X 10 -4 (solid), 2 X 10" 1 (long- dashed), and 1CT 4 (short- dashed). The torque is normalized to 10~ 3 GM s 2 (Mi/M s )2/ai and S is in 
units of M s r^ 2 . Bottom: vertical stratification of the mass density at disk azimuth <f>j = <j>s, in units of M s rj~ , for M2 = 10 — 4 M s (left) 
and 2 X 10 -4 M 3 (right). 



closer to 10~ 4 M S for H/r « 0.03), which implies that 
a substantial fraction of Saturn's envelope is acquired in 
the inner regions of the solar nebula. 

This circumstance, however, would be likely at odds 
with the elemental abundances of some species measured 
in Saturn's atmosphere (see Hersant et al. 2008). The 
abundances relative to hydrogen of elements such as C, 
N, S, As, and P are a few to several times as high as 
the solar abundances (see Lodders 2003; Asplund et al. 
2009, and references therein). In fact, the presence in 
large amount of these elements is believed to have arisen 
from accretion of gas (e.g., Guillot & Hueso 2006) and/or 
solids (e.g., Hersant et al. 2008) in a cold disk environ- 
ment. 

8. SUMMARY AND DISCUSSION 

This paper presents results of thermodynamical mod- 
els of protoplanetary disks that are constructed by ap- 
plying ranges of parameters that may have characterized 
the early solar nebula (see Section 2). Disk evolution is 
driven by viscous torques and photo-evaporation origi- 
nating from the central star (see Section 2.1). Thermal 
balance in the disk is achieved by equating viscous and 
stellar irradiation heating with radiative cooling in the 
vertical direction (see Section 2.2). Only models that 
predict disk lifetimes between 1 and 20 Myr are consid- 
ered viable representations of the solar nebula (see Ta- 
ble 1), in line with observations and core- nucleated ac- 
cretion calculations of gas giants. Such models provide 
the physical conditions (see, e.g., Figure 3) at the time 
and after the planets acquired most of their mass (see 
Section 2.4) and can be used to simulate their long-term 
orbital migration. 

Two and three dimensional hydrodynamical calcula- 
tions are used to quantify the migration rates of a pair of 
planets with mass ratios corresponding to Jupiter's and 
Saturn's, M ± /M s « 10" 3 and M 2 /M s « 3 x 10" 4 (see 



Section 3). The orbits of the planets are initially placed 
in proximity of the 3:2 mean motion resonance (see, e.g., 
Figure 7). As described by Masset & Snellgrove (2001), 
a necessary condition to activate outward migration is 
the overlap of the tidal gap carved in the disk by a less 
massive, exterior planet with the gap of the more mas- 
sive, interior planet (see Figure 8). The relative depth 
and width of the gaps provide the sufficient condition 
(see Section 3.2). High temperatures and kinematic vis- 
cosities inhibit gap formation, either promoting inward 
migration (see Section 3.3) or stopping outward migra- 
tion (see Section 3.4) 

If a near 3:2 commensurability is preserved, there are 
stalling radii in the disk, toward which the pair will con- 
verge (see Figure 14). In general, to first approximation, 
these radii depend on a combination of the turbulence 
viscosity parameter, disk thickness, and mass of the outer 
planet (sec Figures 15 and 22). 

For planets moving outward from the inner disk region 
(r < 2 AU) at a time between 1 and 3 Myr, the interior 
planet may reach beyond ~ 5 AU only if viscosity is low 
enough, the surface density is not too steep, and the disk 
not too warm (see Section 4). However, the probability 
for this to happen appears low (see Table 2). Experi- 
ments performed on random samples of the parameter 
space suggest that in 98% of the cases the interior planet 
stops within 4AU (see Figure 18). 

At least three requirements must be satisfied to estab- 
lish and maintain the 3:2 commensurability and hence 
promote outward migration: 1) the exterior planet must 
stop growing (see Section 5), 2) the interior planet must 
do the same (see Section 6), and 3) the relative migra- 
tion speed prior to capture must be large, so that the 
exterior planet can transit the 1:2 orbital resonance (see 
Section 7). If requirement 1) is violated, the mass ra- 
tio MijM\ may approach 1 on a timescale shorter than 
the migration timescale (see Figure 19), and the outward 
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motion is interrupted. If requirement 2) is violated, the 
accretion rate through the disk, past the interior planet, 
is reduced. The surface density inside the orbit of the in- 
ner planet drops and so does the positive Lindblad torque 
exerted on the planet, inhibiting outward migration (see 
Figure 20). If requirement 3) is violated, migration is 
not reversed and both planets continue moving toward 
the star (see Figure 21). The 1:2 orbital resonance may 
still induce outward migration, but at masses larger than 
Jupiter's and Saturn's (see Figure 16). 

The outward migration mechanism is operable, as also 
argued in previous studies, but the limitations can be se- 
vere. In particular, it is difficult to reconcile the absence 
of accretion on both giant planets with the presence of 
gas around the planets (see, e.g., Lissaucr ct al. 2009, 
and references therein) . Two processes capable of shut- 
ting down the accretion of gas must be invoked although, 
in principle, they need not be different. Additionally, 
since envelope collapse begins once the envelope mass 
exceeds the core mass, it is reasonable to assume that 
before growth stops both planets were undergoing run- 
away gas accretion, i.e., digesting all the gas the nebula 
could provide. Presumably, the sought processes are not 
"internal", i.e., related to the structure of the envelopes, 
because otherwise they would likely occur around sim- 
ilar envelope masses 9 , which is obviously not the case. 
But if the processes are of an "external" nature, i.e., re- 
lated to the supply of gas, then the interior planet would 
probably undergo through said process before the exte- 
rior planet would, since disk gas removal within several 
AU proceeds from the inside out. Therefore, it seems as 
though the problem of stopping gas accretion on both 
planets does not admit a trivial solution. 

The transit of the exterior planet across the 1:2 com- 
mensurability with the interior planet, within a few AU 
of the star, requires surface densities far in excess of those 
predicted by the disk evolution models constructed here. 
The 1:2 resonant-orbit configuration is then favored, in 
a statistical sense, over the 2:3 one. It is worth men- 
tioning that core-nucleated accretion models necessitate 
high enough surface densities of solid material, but in the 
form of planctcsimals not of dust. Since it takes time to 
turn dust into planetesimals and the gaseous disk evolves 
during that time, dust-to-gas mass ratios do not provide 
useful information about gas densities at the time giant 
planets acquired most of their gaseous contents. Disk- 
limited accretion rates depend linearly on the gas sur- 
face density, but tend to be rather large. At ~ 5 AU, if 
£ was of order 10 g cm~ 2 around the time of the runaway 
gas accretion phase, it would take ~ 10 5 years to deliver 
over one half of the current mass of Jupiter (see, e.g., Lis- 
sauer et al. 2009, and references therein). At ~ 9AU, a 
few gem -2 would be sufficient to deliver about a Saturn 
mass worth of gas in ~ 10 5 years. Hence, it is reasonable 
to assume that low gas densities in the solar nebula do 
not prevent the giant planets from reaching their final 
masses. 

Even though it appears unlikely that Saturn can tran- 
sit the 1:2 commensur ability with Jupiter in an evolved 

9 The fast contraction phase that initiates runaway gas accre- 
tion is not much influenced by boundary conditions, i.e., by the 
thcrmodynamical state of the disk. Furthermore, given the com- 
pact orbital configuration of the planets, it is unlikely that disk 
conditions would be very different at the two locations. 



nebula, there is the possibility that it forms within this 
orbital resonance. In such case, however, we argue (see 
Section 7.1) that it may be difficult for the pair to reach 
the 1-2 AU disk region. In fact, capture in the 2:3 mean 
motion resonance with Jupiter and ensuing migration re- 
versal can occur before Saturn attains its full mass (see 
Figure 22), probably when it has between about 1/3 and 
2/3 of the final mass. 



We thank the referee for a prompt response and 
helpful suggestions. G.D. thanks Los Alamos Na- 
tional Laboratory for its hospitality. G.D. acknowl- 
edges support from NASA Outer Planets Research Pro- 
gram grant 202844.02.02.01.75 and from NASA Origins 
of solar systems Program grants NNX11AD20G and 
NNX11AK54G. Resources supporting this work were 
provided by the NASA High-End Computing (HEC) 
Program through the NASA Advanced Supercomputing 
(NAS) Division at Ames Research Center. 



REFERENCES 



Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2005, MNRAS, 
358, 283 4 

Alibcrt, Y., Mousis, O., Mordasini, C, & Benz, W. 2005, 

Astrophys. J., 626, L57 15 
Armitage, P. J. 2011, ARA&A, 49, 195 2 
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, 

ARA&A, 47, 481 22 
Bath, G. T., & Pringle, J. E. 1981, MNRAS, 194, 967 4 
Bcnvcnuto. O. G., Fortier, A., & Brunini, A. 2009, Icarus, 204, 

752 15 

Beutlcr, G. 2005, Methods of celestial mechanics. Vol. I: Physical, 
mathematical, and numerical principles, ed. Beutler, G. 
(Berlin: Springer.) 13 

Brent, R. P. 1973, Algorithms for Minimization Without 
Derivatives (Prentice-Hall, Inc.) 4 

Bromley, B. C, & Kenyon, S. J. 2011, ApJ, 731. 101 15 

Chambers, J. E. 2009, ApJ, 705, 1206 3 

Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and 

Planetary Sciences, 38, 493 4 
Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 3 
Clarke, C. 2011, The Dispersal of Disks around Young Stars, ed. 

P. J. V. Garcia (The University of Chicago Press), 355-418 2 
Crida, A., Masset, F., & Morbidclli, A. 2009, ApJ, 705, L148 15 
D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 

411 3 

D'Angelo, C, Bate, M. R., & Lubow, S. H. 2005, Mon. Not. R. 

Astron. Soc, 358, 316 6, 7, 9 
D'Angelo, G., Henning, T., & Kley, W. 2003, Astrophys. J., 599, 

548 3, 7 

D'Angelo, G., & Lubow, S. H. 2008, ApJ, 685, 560 (DL08) 10, 

12, 16, 19, 20, 21, 23 
D'Angelo, G., Lubow, S. H., & Bate, M. R. 2006. ApJ, 652, 1698 

10, 12, 18 
Davis, S. S. 2005, ApJ, 627, L153 4 

de Val-Borro, M., Edgar, R. G., Artymowicz, P., Ciecielag, P., 
Cresswell, P., D'Angelo, G., Dclgado-Donatc, E. J., Dirkscn, 
G., Fromang, S., Gawryszczak, A., Klahr, H., Kley, W., Lyra, 
W., Masset, F., Mellema, G., Nelson, R. P., Paardekooper, S., 
Peplinski, A., Pierens, A., Plewa, T., Rice, K., Schafer, C, & 
Speith, R. 2006, MNRAS, 370, 529 7 

Dodson-Robinson, S. E., Bodenhcimer, P., Laughlin, G., Willacy, 
K., Turner, N. J., & Beichman, C. A. 2008, ApJ, 688, L99 15 

Dullcmond, C. P., Hollenbach, D., Kamp, I., & D'Alessio, P. 
2007, Protostars and Planets V, 555 2 

Ferguson, J. W., Alexander. D. R„, Allard, F., Barman, T., 
Bodnarik, J. G., Hauschildt, P. H., Heffner-Wong, A., & 
Tamanai, A. 2005, ApJ, 623, 585 3 

Fernandez, J. A., & Ip, W.-H. 1984, Icarus, 58, 109 1 

Franklin, F. A., Lewis, N. K., Soper, P. R., k, Holman, M. J. 
2004, AJ, 128, 1391 1 

Gammie, C. F. 2001, ApJ, 553, 174 7 

Gladman, B. 1993, Icarus, 106, 247 9, 15 

Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, 

Nature, 435, 466 1 
Gorti, U., Dullcmond, C. P., & Hollenbach, D. 2009, ApJ, 705, 

1237 2 

Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539 2 



24 



D'Angelo & Marzari 



Guillot, T., & Hueso, R. 2006, MNRAS, 367, L47 22 
Hahn, J. M., & Malhotra, R. 1999, AJ, 117, 3041 1 
Haircr, E., N0rsctt, S. P., & Wanner, G. 1993, Solving Ordinary 
Differential Equations I: NonstifT Problems (Springer Series in 
Computational Mathematics. Vol. 8, 2nd ed.) 4, 9 
Haisch, K. E., Lada, E. A., & Lada, C. J. 2001, Astrophys. J., 
553, L153 2 

Hartmann, L., Zhu, Z., & Calvet, N. 2011, ArXiv:1106.3343 3 
Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 
70, 35 4 

Hcrsant, F., Gautier. D., Tobie, G., & Lunine, J. I. 2008, 

Planet. Space Set, 56, 1103 22 
Hubeny, I. 1990, ApJ, 351, 632 3 

Hubickyj, O.. Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 
179, 415 15 

Hueso, R.. & Guillot, T. 2005, A&A, 442, 703 3 

Ivanov, P. B., Papaloizou, J. C. B., & Polnarev, A. G. 1999, 

MNRAS, 307, 79 21 
Kenyon, S. J., & Bromley, B. C. 2009, ApJ, 690, L140 15 
Kley, W. 2003, Celestial Mechanics and Dynamical Astronomy, 

87, 85 19 

Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 7, 8 
Kley, W., Peitz, J., & Bryden, G. 2004, A&A, 414, 735 19 
Lin, D. N. C, & Papaloizou, J. 1986a, ApJ, 307, 395 12 
— . 1986b, Astrophys. J., 309, 846 11 

Lissauer, J. J., Hubickyj, O., D'Angelo, G., & Bodenheimer, P. 

2009, Icarus, 199, 338 15, 16, 23 
Lodders. K. 2003, ApJ, 591, 1220 22 

Lubow, S. H., & D'Angelo, G. 2006, Astrophys. J., 641, 526 

(LD06) 18, 19, 20, 24 
Lubow, S. H., & Ida, S. 2010, Planet Migration, ed. Seager, S. 

(Tucson, University of Arizona Press), 347-371 11 
Lynden-Bcll. D., & Pringle, J. E. 1974, Mon. Not. R. Astron. 

Soc, 168, 603 2, 4, 5 
Malhotra, R. 1993, Nature, 365, 819 1 
— . 1995. AJ, 110, 420 1 
Masset, F. 2000, A&AS, 141, 165 7, 8 

Massct, F., & Snellgrove, M. 2001, MNRAS, 320, L55 9, 19, 20, 
22 

Masset, F. S. 2002, A&A, 387, 605 6 

Menou, K., & Goodman, J. 2004, ApJ, 606, 520 3 

Michtchenko. T. A.. Beauge, C, & Ferraz-Mcllo. S. 2008, 

MNRAS, 391, 215 19 
Mihalas, D., & Weibel Mihalas, B. 1999, Foundations of radiation 

hydrodynamics (New York: Dover, 1999) 3 
Minton, D. A., & Malhotra, R. 2009, Nature, 457, 1109 1 
Morbidelli, A., & Crida, A. 2007, Icarus, 191, 158 11, 16 



Morbidclli, A., Lcvison, H. F., Tsiganis, K., & Gomes, R. 2005, 

Nature, 435, 462 1 
Mordasini, C, Alibert, Y., Klahr, H., & Benz, W. 2011, Detection 

and Dynamics of Transiting Exoplanets, St. Michel 

l'Observatoire, France, Edited by F. Bouchy; R,. Diaz; 

C. Moutou; EPJ Web of Conferences, Volume 11, id.04001. 11, 

4001 15 

Movshovitz. N., Bodenheimer, P., Podolak, M., & Lissauer. J. J. 

2010, Icarus, 209, 616 15 
Miiller, T. W. A.. Kley, W., & Meru, F. 2012, A&A, 541, A123 6 
Mustill, A. J., & Wyatt. M. C. 2011, MNRAS, 413, 554 19 
Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, 

Mon. Not. R. Astron. Soc, 318, 18 6 
Papaloizou, J., & Lin, D. N. C. 1984, ApJ, 285, 818 11 
Pierens, A., & Nelson, R. P. 2008, A&A, 482, 333 19, 20 
Pierens, A., & Raymond, S. N. 2011, A&A, 533, A131 1, 12, 14 
Pollack, J. B., Hollenbach, D., Bcckwith, S., Simonelli, D. P., 

Roush, T., & Fong, W. 1994, ApJ, 421, 615 3, 4 
Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., 

Podolak, M., & Greenzweig, Y. 1996, Icarus, 124, 62 15 
Press, W. H., Teukolsky, S. A., Vcttcrling, W. T., & Flannery, 

B. P. 1992, Numerical recipes in FORTRAN. The art of 

scientific computing (Cambridge: University Press, 2nd ed.) 4 
Pringle, J. E. 1981, ARA&A, 19, 137 3 
— . 1991, MNRAS, 248, 754 18 
Quillcn, A. C. 2006, MNRAS, 365, 1367 19 
Rafikov, R. R., & Do Colle, F. 2006, ApJ, 646, 275 3 
Roberge, A., & Kamp, I. 2010, Protoplanctary and Debris Disks, 

ed. Seager. S. (Tucson, University of Arizona Press), 269-295 2 
Ruden, S. P., & Pollack, J. B. 1991, ApJ, 375, 740 2 
Scott, E. R. D. 2006, Icarus, 185, 72 16 
Shakura, N. I., & Syunyaev, R. A. 1973, A&A, 24, 337 4 
Siess, L., Dufour, E., & Forcstini, M. 2000, A&A, 358, 593 4 
Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 7, 8 
Syer, D., & Clarke, C. J. 1995, MNRAS, 277, 758 21 
Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005. 

Nature, 435, 459 1 
Walsh, K. J., Morbidelli, A., Raymond, S. N., O'Brien, D. P., & 

Mandell, A. M. 2011, Nature, 475, 206 1 
Ward, W. R. 1986, Icarus, 67, 164 12 

Ward, W. R., & Hahn, J. M. 2000, Protostars and Planets IV, 
1135 11 

Weidenschilling, S. J. 1977, Ap&SS, 51, 153 4 
Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 2 
Zhang, H., & Zhou, J.-L. 2010, ApJ, 714, 532 9, 19 



